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Abstract 


Despite  the  advances  in  the  theoretical  analysis,  the 
science  of  aerodynamics  still  relies  heavily  on  the  use  of 
wind  tunnels  to  provide  design  and  performance  data  from 
test  models.  Such  data  must  be  corrected  for  the  influence 
of  the  tunnel  walls.  For  tests  where  the  model  is  small  in 
relation  to  the  tunnel  test  section  area  the  correction 
procedures  are  well-established  and  accurate.  Modern 
aerodynamic  research  is  often  concerned,  however,  with  test 
configurations  where  the  model  is  quite  large  compared  to 
the  tunnel  area,  or  develops  very  high  lift.  The  latter 
category  includes  multi-component  sections.  For  these 
geometries  the  standard  correction  forms  become  increasingly 
inadequate,  and  there  are  many  difficulties  with  the 
interpretation  of  such  experimentally  obtained  data. 

A  potential  flow  method  for  determining  the  wind  tunnel 
corrections  has  been  developed  which  contains  explicitly  the 
effect  of  the  wall  interference.  Consideration  was  limited 
to  the  case  of  two-dimensional,  incompressible  flow  in 
closed  wind  tunnels.  There  are  no  restrictions  on  model  size 
or  tunnel  geometry.  The  method  relies  on  a  numerical 
solution  for  the  stream  function  describing  the  complete 
flow  field  in  the  wind  tunnel.  The  solution  is  obtained 
using  finite  difference  techniques  in  a  boundary  fitted 
coordinate  system.  The  coordinate  system  is  generated 
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numerically,  and  allows  an  accurate  and  efficient  solution 
for  the  potential  flow.  Several  types  of  coordinate  systems 
were  investigated  to  determine  the  applicability  of  each  to 
the  wind  tunnel  problem.  With  suitable  grids,  the  numerical 
solution  for  the  entire  flow  field  yields  results  which  are 
comparable  in  accuracy  to  other  potential  flow  analysis 
techniques . 

The  corrections  are  obtained  by  comparison  of  the 
aerodynamic  force  coefficients  obtained  using  the  numerical 
analysis  with  values  obtained  using  a  proven  singularity 
method  for  the  free  flow  field,  that  is,  without  the 
interference  of  the  walls.  The  corrections  so  obtained 
represent  the  combined  effect  of  lift  interference  and  solid 
blockage,  the  coupling  of  which  is  essentially  ignored  by 
the  standard  correction  forms. 

Three  airfoil  sections  of  markedly  different  geometries 
and  performance,  including  one  two-element  section,  were 
tested  in  the  University  of  Alberta  low-speed  wind  tunnel. 
The  tunnel  size  was  made  adjustable  by  using  two  large  wall 
inserts,  allowing  a  wide  range  of  test  geometries  for  each 
model.  The  numerically  generated  corrections  were  found  to 
be  accurate  in  reducing  the  measured  data  to  a  common  set  of 
data  indicative  of  free  field  performance. 

The  procedure  for  obtaining  the  corrections  is  both 
expensive  and  time  consuming,  and  requires  major  computing 
facilities.  Therefore,  no  attempt  was  made  to  develop 
universal  correction  formulae,  since  the  corrections  appear 


v 


to  be  highly  model  dependent  for  these  extreme  geometries. 
However,  the  corrections  are  shown  to  be  capable  of 
substantially  increasing  the  range  of  allowable  test 
configurations  for  which  reliable  wall  corrections  can  be 
made . 
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I.  Introduction 


The  science  of  aerodynamics,  and  its  utilization  in 
practical  aviation,  has  advanced  rapidly  since  the 
introduction  of  the  boundary  layer  theory  in  1904.  This 
major  step  in  the  theoretical  analysis  allowed  the  division 
of  the  fluid  flow  field  into  two  distinct  regions:  the  thin 
boundary  layer  near  solid  surfaces  where  viscous  effects 
could  not  be  neglected,  and  the  external  inviscid  flow  which 
could  be  analyzed  using  the  established  procedures  of 
classical' hydrodynamics .  However,  the  theoretical  analysis 
of  real  fluid  flows  still  could  not  accurately  provide  all 
the  necessary  data  for  aerodynamic  design  purposes,  and 
knowledge  of  the  real  flow  behaviour  could  be  gained  only  by 
experimental  observation.  To  a  large  extent,  the  need  for 
experimental  studies  still  remains,  although  the  development 
of  the  various  mathematical  techniques  has  allowed  a  large 
reduction  in  the  required  experimental  effort.  In  the  last 
three  decades,  an  explosive  growth  in  the  use  of  the 
computer  for  solving  flow  problems  has  occurred. 

Using  mathematical  procedures  developed  much  earlier 
but  impractical  at  the  time  and  by  devising  new  methods, 
theoreticians  can  now  model  many  real  fluid  flows  with 
reasonable  accuracy.  The  theory  and  application  of  numerical 
analysis  techniques  for  problems  in  fluid  flow  has  come  to 
be  called  computational  fluid  dynamics.  Computational 
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methods  presently  form  a  large  portion  of  basic  research  in 
aerodynamics.  Coupled  with  the  ever-increasing  capability  of 
modern  digital  machines,  this  effort  has  both  improved  the 
accuracy  and  increased  the  efficiency  of  numerical 
solutions,  and  allowed  application  to  a  wider  range  of  flow 
problems.  Such  numerical  investigations,  in  what  have*  been 
termed  electronic  wind  tunnels,  can  provide  quick  solutions 
to  analysis  problems  in  which  dozens  of  parameters  are 
involved,  making  it  possible  to  quickly  modify  and  optimize 
a  design  to  meet  particular  goals. 

However,  numerical  methods  are  subject  to  sometimes 
very  severe  restrictions,  either  in  a  mathematical  sense  or 
a  practical  one.  The  physics  of  fluid  motion  is  very 
complex,  and  there  are  no  adequate  mathematical  models  which 
will  handle  the  general  problems  of  aerodynamics.  Often 
these  numerical  simulations  can  be  accomplished  only  by 
making  severe  assumptions  about  the  flow  field.  There  are 
many  physical  phenomena  (turbulent  flow,  separated,  and 
three  dimensional  flows,  to  name  a  few)  for  which  universal 
and  adequate  models  do  not  yet  exist  or  cannot  be 
practically  implemented  due  to  the  limited  capacity  of 
digital  hardware.  It  is  often  these  phenomena  which  are  the 
most  important  features  of  the  flow  and  until  totally 
adequate • numerical  models  and  techniques  are  developed  the 
need  for  experimental  investigation  will  persist. 

The  Wright  brothers  relied  on  aerodynamic  data  obtained 
in  their  own  wind  tunnel,  and  there  are  few  aircraft 
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companies  today  which  will  proceed  with  a  major  project 
without  very  careful  experimental  verification  of  designs. 
The  major  tool  of  these  investigations  is  the  wind  tunnel, 
which  allows  precise  control  over  the  environment  in  which 
an  aerodynamic  model  is  tested.  The  validity  of  the  wind 
tunnel  tests  is  assured,  i.e.  the  behaviour  of  the  full 
scale  design  may  be  extrapolated  with  conf idence 'f rom  model 
tests,  if  certain  conditions  are  satisfied.  When  the  model 
shape  is  scaled  accurately,  these  requirements  are  met  by 
matching  the  Reynolds  number  and  Mach  number  of  the  test 
with  the  full  scale  situation.  This  is  the  well  known 
concept  of  dynamic  similarity. 

Besides  the  attainment  of  dynamic  similarity,  the  wind 
tunnel  and  its  operation  must  be  thoroughly  understood  if 
the  data  it  provides  are  to  be  interpreted  correctly.  This 
is  complicated  by  the  wide  variety  of  tunnel  configurations 
and  sizes,  and  the  diverse  nature  of  the  experimental  tests 
which  might  be  performed.  Furthermore,  the  wind  tunnel  does 
not  provide  a  free  flow  field  in  which  the  test  will  be 
carried  out.  Only  rarely  are  tests  carried  out  for  bounded 
flows,  such  as  testing  for  ground  effect  of  aircraft  or  air 
cushion  vehicles.  The  fact  that  the  wind  tunnel  flow  field 
is  bounded  can  affect  the  data  as  much  as  variations  in 
Reynolds  or  Mach  number.  These  boundary  constraints  on  the 
flow  must  be  clearly  understood,  and  their  effects  carefully 
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There  are  two  approaches  to  the  problem  of  boundary 
constraint  or  tunnel  wall  interference.  Usually,  the  effect 
of  the  boundaries  is  compensated  by  making  corrections  to 
the  measured  data.  These  standardized  corrections  are  found 
either  empirically  or  using  simple  mathematical  theories, 
based  on  several  assumptions  about  the  tunnel  geometry.  The 
second  approach,  the  subject  of  much  recent  research,  is  the 
design  of  tunnels  and  test  geometries  which  are  interference 
free.  This  may  be  done  by  varying  the  shape  of  the  tunnel 
walls  or  providing  semi-open  walls,  so  that  the  interference 
is  minimized  or  totally  eliminated.  Such  modifications  would 
need  to  be  based  on  very  large  and  expensive  research 
efforts  in  their  own  right.  The  prospects  for  modification 
to  an  interference-free  configuration  seem  minimal,  assuming 
such  modifications  could  be  physically  accommodated,  if  the 
tunnel  is  to  host  a  wide  range  of  test  geometries. 

Much  recent  research  has  been  directed  towards  high 
lift  systems  operating  at  low  Mach  numbers  and  high  Reynolds 
numbers,  such  as  occur  during  take-off  and  landing  of  large 
aircraft.  To  obtain  the  large  Reynolds  numbers,  ever  larger 
models  are  being  tested  which  press  to  the  limit  the 
applicability  of  the  standard  correction  techniques.  Very 
little  high  speed  work  is  done  with  large  models.  For  these 
tests,  the  Mach  number  is  of  primary  importance,  with  the 
Reynolds  number  assuming  a  secondary  role.  Because  of  power 
demands,  tests  where  compressibility  is  important  are 
therefore  usually  carried  out  using  small  models  in  special 
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The  aim  of  the  present  work  is  to  investigate  the  use 
of  certain  numerical  techniques  in  extending  the  range  of 
applicability  of  the  standard  method  of  applying  boundary 
corrections  for  non-standard  geometries  and  high  lift 
sections.  Although  in  principle  the  numerical  technique 
developed  could  be  extended  to  three  dimensional  flows  and 
open  tunnels,  consideration  will  be  limited  to  the  case  of 
two  dimensional  flow  in  closed  jet  tunnels  with  fixed  walls. 
There  are  several  motivations  for  this.  First,  the  vast 
majority  of  wind  tunnels  operating  today  are  of  this  type, 
and  high  lift  research  is  primarily  concerned  with  the  flow 
past  the  two  dimensional  wing  section  profile.  Second,  the 
existing  methods  for  obtaining  wall  corrections  are  valid 
only  for  a  limited  range  of  test  geometries  which  will  be 
detailed  in  the  next  chapter.  Finally,  the  investigation  of 
new  and  advanced  numerical  procedures  in  aerodynamic 
research  will  help  to  determine  the  applicability  of  such 
methods.  Wind  tunnel  testing  can  become  very  expensive  and 
time  consuming  and  although  computational  methods  will  not 
likely  replace  the  wind  tunnel,  they  will  undoubtedly 
continue  to  assume  a  larger  share  of  the  work  load.  It  is 
worthwhile  therefore  to  explore  the  use  of  these  numerical 
techniques,  to  identify  useful  applications  and  areas  for 
improvement . 

As  pointed  out  by  Garner  in  [1],  there  is  a  need  for 
improved  understanding  of  the  wall  interference  effects  when 
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the  model  is  large  in  relation  to  the  tunnel  cross  section 
and/or  develops  high  lift.  In  these  cases,  the  problem 
becomes  non-linear  and  the  various  wall  effects  interact. 
Garner  suggests  that  ([1],  p.15)  "a  convincing  theory  of 
wall  interference  on  high  lift  systems  may  require  an 
improved  mathematical  model"  and  that  ([1],  p.14)  "the  most 
promising  method  of  solution  may  well  be  one  of  numerical 
analysis  by  a  finite-difference  technique."  The  coupled 
effects  of  the  tunnel  walls  may  then  be  calculated 
simultaneously . 

The  present  work  is  an  investigation  of  just  such  a 
solution  scheme.  By  eliminating  many  of  the  underlying 
assumptions  used  to  develop  the  standard  correction  schemes, 
the  analysis  results  in  corrections  which  are  applicable  to 
a  much  wider  range  of  test  situations.  In  the  case  of 
multi-element  and  high  lift  sections,  the  present  analysis 
offers  an  accurate  and  intuitively  correct  interpretation  of 
the  wall  interference  corrections. 

The  theoretical  analysis  is  accomplished  using  a  finite 
difference  solution  for  the  stream  function  which  describes 
the  tunnel  flow.  The  solution  for  a  given  test  geometry  is 
carried  out  in  a  numerically  generated  boundary  fitted 
coordinate  system.  This  coordinate  system  offers  the 
advantage  of  placing  computational  nodes  exactly  on  the  flow 
boundaries  and  distributing  the  nodes  in  a  rational  manner 
throughout  the  field.  This  technique  alleviates  the 
difficulty  of  performing  finite  difference  solutions  on 
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fields  with  irregular  boundaries.  For  the  present  work, 
several  grid  topologies  were  investigated,  to  determine  the 
applicability  of  each  grid  type  to  the  given  flow  problem. 
Each  grid  topology  offered  certain  advantages  to  the  present 
solutions.  By  carefully  controlling  the  nature  of  the 
computational  grid  in  the  physical  plane,  adequate 
theoretical  results  were  attainable.  These  results  lead 
directly  to  the  estimation  of  the  corrections  necessary  to 
account  for  the  interference  of  the  walls  during  the  wind 
tunnel  tests. 

In  the  following  chapters  the  material  is  organized  as 
follows.  Chapter  II  reviews  the  existing  theories  for  wind 
tunnel  boundary  corrections.  The  need  for  a  more  general 
theory  will  be  pointed  out.  The  following  chapter  describes 
the  theoretical  approach  developed  to  handle  arbitrary  test 
configurations.  By  comparing  solutions  for  the  bounded  wind 
tunnel  flow  with  those  for  a  free  field  flow,  the 
interference  of  the  walls  may  be  assessed,  and  suitable 
corrections  devised  for  each  tunnel  geometry.  Chapter  IV 
describes  some  experiments  carried  out  in  the  University  of 
Alberta  low  speed  wind  tunnel  to  demonstrate  the 
applicability  and  accuracy  of  the  proposed  corrections,  with 
a  comparison  with  standard  forms,  for  a  wide  range  of  test 
conf igurat ions . 

The  results  indicate  that  the  present  method  can 
generally  provide  better  corrections  than  the  standard 
methods  for  the  extreme  geometries  tested.  Thus,  the 
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permissible  range  of  geometries  for  which  reliable 
corrections  are  available  has  been  extended.  The  method 
applies  equally  well  to  multi-element  and  high  lift 
sections,  which  are  not  adequately  treated  by  the  standard 
corrections . 


II.  Wind  Tunnel  Corrections  -  Standard  Forms 


Experimental  investigations  in  wind  tunnels  are  usually 
directed  towards  the  measurement  of  the  net  forces  and 
moments  exerted  on  the  model  by  the  moving  stream.  This  data 
is  to  be  applied  to  the  estimation  of  the  full  scale 
performance  which  occurs  in  the  free  air  of  the  atmosphere. 
There  are  two  fundamental  concerns  associated  with  the 
generation  of  reliable  wind  tunnel  test  data.  The  first  is 
the  attainment  of  dynamic  similarity,  i.e.  matching  the 
Reynolds  and  Mach  numbers,  which  ensures  that  the  wind 
tunnel  flow  corresponds  geometrically  and  dynamically  with 
the  full  scale  flow.  The  second  concerns  the  estimation  of 
the  effect  of  the  tunnel  boundaries  on  the  flow  surrounding 
the  model,  specifically  the  effect  on  the  measured 
aerodynamic  forces  and  moments.  The  nature  of  these  effects 
depends  upon  whether  the  tunnel  is  of  the  open  jet  type,  or 
constrains  the  flow  between  rigid  walls.  (There  are  numerous 
other  experimental  considerations,  such  as  free  stream 
turbulence  and  non-uniformity,  growth  of  wall  boundary 
layers,  and  the  presence  of  measuring  gear  in  the  test 
section,  which  will  not  be  examined  here.) 

The  boundary  effects  are  usually  obtained  as  small 
incremental  corrections  which  must  be  applied  to  the 
measured  data  such  that  the  corrected  data  is  indicative  of 
the  free  field  performance.  These  incremental  adjustments 
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depend  upon  the  geometry  of  the  wind  tunnel  and  model, 
primarily  the  chord-to-height  ratio.  For  usual  geometries 
the  corrections  may  be  5%  to  10%  of  the  measured  values. 
Various  theoretical  approaches,  which  will  be  discussed 
shortly,  have  been  developed  for  the  estimation  of  these 
corrections . 

'  The  methods  to  be  treated  here  are  concerned  with  two 
dimensional  testing  in  tunnels  with  closed  walls.  Very 
little  two-dimensional  testing  is  done  in  open  tunnels 
because  the  flow  is  often  far  from  two-dimensional  due  to 
end  effects.  However,  the  open-wall  corrections  are  of  the 
opposite  sign  to  the  closed-wall  corrections,  which  raises 
the  possibility  of  providing  semi-open  walls  where  the 
corrections  should  cancel.  Such  a  configuration  was 
developed  by  Williams  and  Parkinson  [2],  who  obtained  a  zero 
correction  tunnel  using  a  slotted  upper  wall  with  an  open 
area  ratio  of  0.7.  With  the  slotted  wall  extending  from  two 
to  six  chord  lengths  upstream  and  downstream  from  the  model, 
application  of  this  method  seems  limited  for  large  tunnels 
due  to  physical  constraints. 

Other  zero  interference  tunnels  have  been  described  by 
Goodyer  [3]  and  by  Vidal,  Erickson,  and  Catlin  [4],  The 
self-streamlining  tunnel  of  Goodyer  utilizes  flexible  upper 
and  lower  walls,  which  are  automatically  positioned  during 
the  test  to  simulate  streamlines  in  a  free  flow  field.  This 
technique  does  not  appear  to  be  applicable  to  large  tunnels, 
again  because  of  the  extent  of  the  required  modifications. 
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The  self-correct ing  tunnel  of  Vidal,  et  al.,  utilizes  large 
regions  of  suction  and/or  blowing  through  porous  walls,  the 
result  being  the  duplication  of  the  free  field  streamline 
pattern  at  the  model.  This  method  was  developed  for 
transonic  flow,  where  interference  of  the  wall  with  the 
shock  waves  is  critical,  and  may  be  impractical  for  large 
tunnels  due  to  power  demands. 

Figure  1  represents  diagramat ically  the  situation  which 
occurs  in  a  closed  tunnel  with  straight  walls,  where  the 
measured  lift  coefficient  C^.  at  the  test  angle  of  attack  0L,. 
is  compared  to  the  free  field  data.  It  is  seen  that  the 
measured  lift  is  too  high,  the  measured  incidence  is  too 
low,  or  both.  The  measured  drag  and  pitching  moment  are  also 
larger  than  free  field  values.  The  three  lines  I,  II,  and 
III,  represent  three  possible  correction  strategies  to 
obtain  proper  free  air  lift  coefficients,  C^_  at  9p  ,  Cjp  at 
9p  ,  etc.  These  strategies  will  be  discussed  at  length  in  a 
later  section.  For  the  present,  it  is  necessary  to  identify 
the  sources  of  the  experimental  deviations. 

Theoretically,  the  need  for  the  tunnel  wall  corrections 
may  be  deduced  from  the  fact  that  although  the  governing 
flow  equations  are  the  same  for  the  flow  in  the  tunnel  and 
in  free  air,  the  outer  boundary  conditions  on  the  flow  field 
are  different.  In  free  air,  the  streamlines  far  from  the 
airfoil  are  able  to  curve  freely,  whereas  in  the  tunnel,  the 
flow  is  constrained  between  straight  rigid  walls.  For 
subsonic,  steady  flow  the  governing  equations  are  elliptic 
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Figure  1  Effect  of  Tunnel  Walls  on  Measured  Lift 
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and  the  boundary  conditions  have  a  major  effect  on  the  flow 
solution,  and  hence,  on  the  aerodynamic  forces  on  the  model. 

The  increments  in  the  force  coefficients  caused  by  wall 
interference,  hence  the  corrections  required  to  the 
experimental  data,  may  be  attributed  to  two  major  effects. 
These  are  described  in  Pope  and  Harper  -[5]  and,  in  greater 
detail,  in  Garner,  Rogers,  Acum  and  Maskell  [1].  Of  primary 
importance  to  two  dimensional  testing  in  closed  tunnels  are 
the  corrections  required  to  account  for  the  interference  due 
to  the  blockage  of  the  tunnel  by  the  model  and  its  wake,  and 
to  the  constraining  effect  of  the  straight  walls  on  the 
curvature  of  the  flow  around  the  model.  The  geometric 
parameters  which  govern  these  corrections  are  defined  in 
Figure  2.  These  may  be  grouped  into  the  non-dimensional 
terms  representing  the  chord-to-height  ratio  (c/h),  model 
position  (d/h) ,  and  model  rotation  point 
(  Xr  /c).  Along  with  the  angle  of  attack  (9r(and  the  flap 
angle  for  two  element  sections),  these  determine  the 
basic  tunnel  geometry.  The  profile  of  the  model  itself  may 
be  described  in  rough  terms  by  the  thickness-to-chord  ratio 
(t/c),  cross-sectional  area  (A),  and  the  camber  ratio  (y/c). 

Early  work  on  the  theory  of  tunnel  corrections,  being 
analytical  in  nature,  was  possible  only  by  restricting  the 
consideration  to  the  case  of  small,  thin  models  at  low  angle 
of  attack  whose  mid-chord  was  on  the  tunnel  centre-line.  The 
tunnel  was  also  assumed  to  extend  infinitely  far  upstream 
and  downstream  from  the  model.  Later  theories  included 
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Figure  2  Definition  of  Tunnel  Geometry 
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special  cases  for  off-center  model  locations,  and 
consideration  of  the  effects  of  thickness  and  larger  c/h 
ratios . 

For  incompressible  unseparated  flow,  the  wall  effect 
may  be  analysed  using  potential  flow  techniques,  by  ignoring 
the  effects  of  viscosity.  Historically,  the  most  important 
of  these  is  the  method  of  images,  Figure  3,  where  the 
influence  of  the  walls  is  calculated  from  the  flow  field 
induced  at  the  model  by  the  infinite  set  of  mirror  images 
above  and  below  the  tunnel  (only  the  first  images  are 
shown).  Of  necessity  these  early  analyses  were  limited  to 
single  airfoils,  and  mathematically  soluble  only  by  making 
the  assumptions  described  previously. 

Extensions  to  these  early  theories  'have  been  developed 
to  account  for  higher  order  corrections  for  larger  models, 
as  described  below.  These  methods  have  been  used  for  test 
cases  at  moderate  chord-to-height  ratios,  but  due  to 
numerical  complexity  are  not  in  common  use. 

In  arriving  at  corrections  using  the  method  of  images, 
the  implicit  assumption  was  made  that  the  tunnel  blockage 
and  streamline  curvature  effects  were  independent.  Thus, 
since  the  governing  flow  equation  is  linear,  the  blockage  : 
and  curvature  effects  could  be  analysed  separately,  using 
image  systems  composed  of  either  doublets  or  vortices.  More 
advanced  image  theories  used  continuous  source-sink  or 
vorticity  distributions.  In  either  case,  the  combined  wall 
effect  was  obtained  by  superposition.  Following  this 
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Figure  3  Method  of  Images 
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procedure,  the  wall  effects  will  be  discussed  separately, 
although  it  develops  (section  D)  that  the  corrections  are 
not  independent  except  under  special  circumstances.  In  the 
discussion  which  follows,  the  subscripts  "T"  and  "F"  refer 
to  measured  test  values  and  corrected  free  field  values, 
respectively . 


A.  Tunnel  Blockage 

Tunnel  blockage  corrections  arise  from  the  volume 
occupied  by  the  model  and  its  wake  in  the  tunnel.  The  model 
reduces  the  cross-sect ional  area  of  the  tunnel,  and  since 
the  flow  is  incompressible,  the  effective  free  stream 
velocity  at  the  model  is  increased  slightly  above  the  tunnel 
speed  measured  far  upstream.  The  required  correction  to  the 
tunnel  velocity  is  usually  divided  into  independent  solid 
and  wake  blockage  components,  where  the  former  refers  to  the 
physical  size  of  the  model  and  the  latter  relates  to  the 
displacement  thickness  of  the  airfoil  boundary  layers  and 
wake,  and  hence  is  related  to  the  airfoil  drag.  The  ratio  of 
the  total  velocity  increment  to  the  free  stream  velocity  is 
represented  by 


(1) 


where  and  £w  are  the  solid  and  wake  blockage  velocity 
increments,  respectively. 
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Solid  blockage 

The  solid  blockage  component  depends  on  the  model 
chord,  maximum  thickness,  and  thickness  distribution,  and  is 
considered  to  be  independent  of  camber  and  lift.  Thus,  the 
determination  of  is  based  on  consideration  of  the  section 
thickness  distribution  obtained  with  the  camber  removed,  and 
mounted  at  zero  incidence.  Early  analyses  were  also  based  on 
the  assumption  that  only  the  longitudinal  flow  was  affected, 
which  is  only  true  if  the  model  is  mounted  on  the  tunnel 
center  line.  This  early  work  is  well  documented  in  Garner, 
et  al.  [1].  Briefly,  these  solid  blockage  corrections  were 
developed  as  follows: 

1.  Lock  (1929)  represented  the  airfoil  and  images  by 

doublets.  The  blockage  factor  obtained  was  given  by  the 
formula 


(2) 


where  X2  r  the  body-shape  factor,  is  a  function  of  the 
thickness  distribution,  values  of  which  are  tabulated 
for  various  sections. 

2.  Goldstein  [6],  using  a  series  solution  technique, 

developed  a  formula  accurate  for  higher  c/h  ratios.  This 
formula  was  later  re-arranged  by  Allen  and  Vincenti  [7] 


to  the  form 
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where  A  is  the  cross  sectional  area  of  the  profile,  and 
(x,y)  are  the  coordinates  of  the  thickness  distribution. 
The  variable  4>  is  given  by  the  function 

x  =  j  (1  -  cos4>) 


3.  Thom  (1943)  replaced  the  doublet  of  Lock  by  a 

distribution  of  sources  and  sinks  on  the  chord  line  to 
obtain  the  formula  (for  thin  sections) 


E  =  =  0.52  \  (4) 

s  6h‘  n 

4.  Young  and  Squire  (1945)  proposed  that  the  formula  of 
Thom  could  be  improved  for  thicker  sections  using 


(5) 


5.  Thompson  (1948)  proposed  a  different  improvement  to 
Thom's  relation  by  including  a  term  which  depended 
explicity  on  the  model  thickness,  making  the  formula 
valid  up  to  t/c=0.2,  as  follows, 


+ 
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More  complex  expressions  were  developed  to  account  for 
streamwise  variations  in  blockage,  Woods  (1954),  and  for 
variations  due  to  angle  of  attack  and  off-center  model 
position,  Batchelor  (1944).  The  latter  correction  is 
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For  small  models  at  low  lift  coefficients  these 
corrections  have  been  proven  to  be  accurate  in  numerous 
experimental  tests.  Allen  and  Vincenti  [7]  showed  the 
accuracy  of  the  Goldstein  corrections  up  to  c/h=0.4  at  low 
Mach  numbers. 

Thus,  for  the  cases  of  large  models,  high  lift 
coefficients,  or  off-center  mounting  position,  the  above 
corrections,  which  become  quite  large,  are  questionable.  In 
these  cases,  such  effects  as  the  chordwise  variation  in  the 
longitudinal  velocity  increment  and  induced  transverse 
velocities  must  be  considered.  These  will  interact  with  the 
corrections  to  lift  and  incidence  described  in  a  later 
section . 


Wake  Blockage 

The  blockage  effect  of  the  wake  behind  the  model  is 
obtained  by  representation  of  the  wake  by  a  source  at  the 
section  trailing  edge  and  a  sink  far  downstream.  The 
velocity  increment  at  the  model  due  to  the  source-sink  image 
system  may  then  be  related  to  the  original  source  strength, 


and  ultimately  to  the  airfoil  drag  coefficient.  Using  this 
approach,  the  wake  blockage  increment 
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has  been  developed  by  several  authors,  for  example  Allen  and 
Vincent!  [7],  Goldstein  [6]  showed  in  addition  that  the 
resulting  correction  is  valid  when  the  drag  is  small  even  if 
the  chord  to  height  ratio  is  large.  The  wake  blockage 
increment  is  much  smaller  than  the  solid  blockage  increment 
as  long  as  the  drag  remains  small,  i.e.  the  flow  remains 
attached  to  the  airfoil  surface.  For  the  present  work,  which 
utilizes  inviscid  theory  and  therefore  restricts 
cons-iderat ion  to  unseparated  flow,  the  above  wake  blockage 
formula  will  be  assumed  to  be  adequate. 

B.  Lift  Interference 

The  presence  of  the  walls  constrains  the  natural 
curvature  of  the  flow  around  a  lift-producing  airfoil.  The 
resulting  corrections  to  the  stream  direction  and  curvature 
are  known  as  lift  interference  corrections  since  they  arise 
solely  from  the  lift  or  circulation  generated  by  the  model. 
The  effect  of  the  interference  is  to  increase  the  effective 
camber  and  incidence  of  the  model,  while  there  is  no  effect 
on  the  longitudinal  velocity  if  the  model  is  mounted  on  the 
tunnel  center-line.  To  a  first  approximation,  the  lift 
interference  is  independent  of  the  model  thickness,  and  so 
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can  be  analysed  independently  of  the  blockage.  However,  it 
is  important  even  for  small  models,  vanishing  only  at  zero 
lift. 

The  methods  of  analysis  detailed  below  ignore  the  model 
thickness  and  place  either  the  mid-chord  or  the  aerodynamic 
centre  on  the  tunnel  center-line.  The  model  is  represented 
by  a  single  vortex  or  continuously  distributed  vorticity 
along  the  chord  line  and  the  effect  of  the  walls  is  deduced 
using  the  method  of  images.  The  primary  effects  of  the  image 
system  are  treated  as  an  induced  upwash  and  an  induced  flow 
curvature.  The  upwash  is  evaluated  at  the  model  mid-chord 
and  treated  as  an  incidence  correction.  The  flow  curvature 
produces  a  change  in  effective  camber  or  loading  and  yields 
the  required  correction  to  the  section  lift  and  pitching 
moment . 

A  complete  discussion  of  lift  interference  appears  in 
Garner,  et  al.  [1],  highlights  of  which  appear  below.  The 
theories  are  classed  as  first  or  second  order,  depending  on 
the  terms  retained  in  the  final  expressions  for  the 
corrections.  The  second  order  theories  are  generally 
preferred  when  the  c/h  ratio  becomes  large.  Briefly,  the 
development  of  the  lift  interference  theories  has  occurred 
as  follows: 

1.  Glauert  (1933)  considered  small,  thin  airfoils  at  small 
incidence.  In  developing  this  first  order  theory,  the 
airfoil  was  represented  by  a  single  vortex  at  the  centre 
of  pressure  giving  the  induced  upwash  angle  at  position 
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x  on  the  chord  line  as 
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This  is  one  of  the  few  theories  which  consider 

off-center  models,  the  effect  of  which  is  shown  by  the 

last  term.  This  arises  from  an  induced  longitudinal 
* 

velocity  component  which  vanishes  when  d=h/2. 

Tomotika  (1938)  used  a  conformal  transformation  method 
to  analyse  a  flat  plate  mounted  on  the  tunnel 
center-line,  obtaining  variations  in  lift  and  moment 
with  angle  of  attack  for  various  chord-to-height  ratios. 
Allen  and  Vincenti  [7]  devised  a  first  order  correction 
for  thin  airfoils  by  replacing  the  model  with  a 
continuous  vorticity  distribution  on  the  chord  line.  The 
upwash  at  the  mid-chord  induced  by  the  image  system  is 
given  by 
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This  agrees  with  the  expression  of  Glauert,  equation 
(9),  evaluated  at  x/c=.5  and  d/h= . 5 .  The  induced 
curvature  results  in  lift  and  moment  corrections 
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With  these  corrections  the  new  lift  and  moment 


24 


coefficients  relate  to  the  unconfined  flow  at  the 
corrected  incidence  9y  +  A0  .  Comparison  of  these 
results  with  Tomotika's  results  showed  large 
discrepancies  for  c/h>0.4,  indicating  the  need  for 
higher  order  corrections  for  these  geometries. 

4.  Havelock  (1938)  considered  a  flat  plate  whose  mid-point 
was  an  arbitrary  distance  from  the  tunnel  floor,  and 
neglected  terms  of  order  greater  than  (c/h)^  to  obtain  a 
second  order  method.  Havelock  showed  that  even  for  small 
c/h,  the  effect  of  incidence  and  higher  terms  in  c/h  can 
be  of  primary  importance  when  the  model  is  off  the 
center-line.  Moreover,  the  downward  displacement  of  the 
model  at  positive  incidence  can  change  the  sign  of  the 
correction  due  to  the  induced  longitudinal  velocity, 
which  can  no  longer  be  ignored. 

5.  Goldstein  [6]  approached  the  general  problem  of  the 
.thick  cambered  section  using  power  series  in  (c/h).  The 

derived  second  order  corrections,  valid  only  when  the 
model  is  on  the  center-line,  are  as  follows: 
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Goldstein  showed  that  for  c/h>0.3  the  terms  in 
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which  reduce  the  magnitude  of  the  correction,  cannot  be 
ignored.  Goldstein  also  solved  the  problem  of  the 
off-center  model  using  power  series,  but  it  suffers  from 
very  slow  convergence. 

6.  Moses  [9]  derived  a  conformal  transformation  technique 
to  calculate  the  velocity  distributions  on  an  arbitrary 
airfoil  located  anywhere  within  the  tunnel,  although 
providing  results  only  for  the  center-line  position. 
Results  were  similar  to  Goldstein's  method,  and  in  the 
opinion  of  Garner,  easier  to  obtain.  Also  according  to 
Garner,  the  Moses  analysis  showed  that  if  terms  in  (c/h) 
are  neglected,  the  corrected  lift  curve  slope  is  too 
low,  while  if  the  thickness  is  neglected,  the  slope  is 
too  high,  by  a  few  per  cent. 

As  with  the  solid  blockage  corrections,  there  have  been 
numerous  experimental  verifications  of  the  above  lift 
interference  corrections.  Knechtel  [8]  showed  the  accuracy 
of  the  first  order  correction  of  Allen  and  Vincenti  for  low 
Mach  numbers  so  long  as  the  c/h  ratio  was  less  than  0.4.  For 
incompressible  flow,  the  second  order  correction  of 
Goldstein  is  preferred  when  the  c/h  ratio  is  large. 

C.  Application  of  the  Corrections 

Assuming  that  the  blockage  and  lift  interference 
corrections  are  independent,  they  may  be  applied  separately 
to  the  measured  wind  tunnel  data.  The  blockage  correction 
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results  in  a  change  in  the  tunnel  velocity  at  the  model. 
Thus,  the  tunnel  speed  is  corrected  by 

Vp  *  VT  (1  +  e) 


where VTis  the  tunnel  speed  measured  far  upstream,  or  at 
least  obtained  without  interference  from  the  model.  The 
dynamic  pressure  is  therefore  corrected  using 


qF  =  qT  0  +  e)2 

The  measured  aerodynamic  force  coefficients  are  corrected 
for  blockage  using  the  corrected  dynamic  pressure: 


and  if  higher  order  terms  in  e  are  neglected,  these 
corrections  become 
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The  correction  due  to  the  lift  interference  involves 
one  of  the  three  strategies  shown  in  Figure  1 .  The  usual 
strategy  is  to  correct  both  the  incidence  and  the  lift 
coefficient,  line  I,  using  either  equations  (10)  or  (11). 
Thus,  with  a  defined  by 
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the  first  order  corrections  of  Allen  and  Vincenti,  equations 
(10)  are  ~ 
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and  the  second  order  corrections  of  Goldstein,  equations 
(11),  become 
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The  total  corrections  are  finally  obtained  by  adding 
the  contributions  due  to  blockage  and  lift  interference. 
Thus,  the  complete  first  order  corrections  become 
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These  are  the  corrections  given  in  Pope  and  Harper  [5],  who 
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unfortunately  do  not  stress  the  limitations  of  their  use. 

If  the  model  is  off  the  center-line,  the  corrected 
incidence  at  the  mid-chord  position  becomes,  using  equation 
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Similar  expressions  may  be  developed  for  the  second  order 
corrections. 

The  two  correction  strategies  indicated  by  lines  II  and 
III  in  Figure  1  represent  corrections  at  constant  incidence 
and  constant  lift,  respectively,  as  given  by  Pankhurst  and 
Holder  [10],  Both  of  these  strategies  require  knowledge  (or 


assumption)  of  the  free  field  lift  curve  slope 
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approximately  equal  to  2n  for  thin  airfoils. 

For  the  correction  at  constant  incidence,  line  II,  the 


lift  correction  is  augmented  by  an  additional  increment 
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corresponding  to  the  induced  upwash  angle.  Thus, 


referring  to  Figure  1, 
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The  correction  at  constant  lift,  line  III,  is  accomplished 
similarly  by  treating  the  lift  increment  due  to  the  camber 
effect  as  an  additional  incidence  correction,  thus, 
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Theoretically,  for  both  these  strategies,  the  pitching 
moment  and  drag  coefficients  should  be  adjusted  accordingly. 
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.  For  attached  flow  however,  these 


derivatives  are  usually  very  small  and  the  corresponding 
additional  increments  are  ignored.  For  either  strategy, 
then,  the  corrected  pitching  moment  and  drag  coefficients 
are 
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The  constant  incidence  strategy  will  be  used  for  the 
application  of  corrections  in  the  experimental  work 
described  in  Chapter  IV  as  this  allows  easier  comparison 
with  the  theoretical  results  of  Chapter  III.  This  affects 
only  the  lift  interference,  as  the  blockage  corrections  are 
assumed,  in  the  above  analyses,  to  be  independent  of  lift 
and  incidence.  These  corrections  result  in  the  following 
expressions  (to  first  order) 
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The  second  order  corrections  take  the  form 
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For  thicker  sections,  and  for  multi-element  sections, 
the  theoretical  inviscid  lift  curve  slope  is  larger  by  a 
small  fraction  than  the  thin  airfoil  value  of  2 rr  (per 
radian).  The  actual  lift  curve  slope  obtained  experimentally 
is  slightly  less  owing  to  viscous  effects,  and  is  therefore 
also  affected  by  the  Reynolds  number.  For  moderately  thick 
single  element  sections  the  effects  of  thickness  and 
viscosity  roughly  balance,  so  that  the  real  lift  curve  slope 
is  near  the  value  2n  .  For  very  thick  and/or  multi-element 
sections,  however,  the  situation  is  not  so  clear. 
Experimental  values  much  higher  than  2rr  have  been  reported: 
Wentz  [11]  measured  a  slope  of  1.29*2rr  for  the  NASA  section 
GA(W)-1  with  a  30%  Fowler  flap  at  20° deflect  ion .  Foster, 
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Irwin,  and  Williams  [12]  measured  1.09*2it  for  a  section 
with  a  40%  slotted  flap  at  30° def lect ion .  The  inviscid  two 
element  analysis  of  Kennedy  [13]  calculates  a  slope  of 
1.25*2rr  for  the  slotted  flap  section  used  in  the  present 
experimental  work.  Unfortunately,  the  real  lift  curve  slope 
for  thick  or  multi-element  sections  cannot  be  known  a 
priori .  For  cases  where  the  experimental  lift  curve  slope 
differs  markedly  from  the  value  2rr  ,  the  constant  incidence 
corrections  must  be  applied  in  several  steps.  The  first  step 
utilizes  the  theoretical  value  in  equation  (19b)  to  obtain 
an  approximation  to  the  corrected  lift  curve.  The  slope  of 
this  curve  is  then  calculated  and  the  corrections  are 
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repeated  with  the  improved  value  for _ F#  This  procedure  is 

seF 

carried  out  iteratively,  and  was  used  in  the  correction  of 
the  experimental  data  described  in  Chapter  IV.  Two 
iterations  were  normally  required. 


D.  Viscous  Effects  and  Non-Linear  Coupling 

Except  for  the  wake  blockage,  the  previous  corrections 
are  all  based  on  inviscid  flow  theory.  So  long  as  the  flow 
remains  attached  this  should  be  satisfactory,  with  one 
exception.  When  the  wall  interference  is  large,  the  airfoil 
pressure  distribution  can  be  altered  considerably.  On 
certain  airfoils  where  the  pressure  distribution  has  a 
critical  effect  on  the  boundary  layers  (laminar  flow 
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airfoils),  this  interference  can  drastically  alter  the  model 
behaviour  far  beyond  the  normal  corrections.  The  ability  to 
predict  such  behaviour  in  order  to  estimate  corrections 
requires  a  complete  and  accurate  theory  for  the  viscous  flow 
field  itself  which,  if  available,  precludes  the  need  for  the 
tunnel  tests  in  the  first  place.  Thus,  caution  must  be 
exercised  when  testing  such  critical  boundary  layer  models, 
and  extreme  tunnel  geometries  should  be  avoided. 

When  the  model  is  thick,  large,  or  mounted  off  the 
tunnel  center-line,  the  blockage  and  lift  interference  are 
no  longer  independent.  Blockage  will  then  induce  transverse 
velocity  components,  thereby  affecting  the  incidence  and 
effective  camber  of  the  section.  Lift  interference  will 
induce  longitudinal  components  which  will  interact  with  the 
solid  blockage.  For  multi-element  sections  the  situation 
becomes  much  more  obscure  due  to  the  extreme  flow  curvature 
behind  high  lift  flaps. 

For  tunnel  geometries  where  the  corrections  are 
coupled,  it  is  reasonable  to  formulate  a  correction  theory 
where  a  single  correction  is  calculated  which  simultaneously 
includes  both  effects.  The  theory  of  the  present  work, 
described  in  the  following  chapter  is  one  such  method. 

E.  Modern  Methods 

It  is  observed  that  the  theory  for  wall  corrections  for 
small  thin  sections  at  low  incidence  was  well  established  by 
the  mid-century.  At  that  time,  the  computational  limit  of 
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existing  methods  was  reached.  The  extremely  involved 
conformal  analysis  method  of  Moses  [9]  is  a  case  in  point. 
Without  the  benefit  of  fast  computation  the  method  is  quite 
unwieldy  and  may  explain  why  the  conformal  methods  and 
similar  exact  analyses  for  correction  estimation  are  not 
more  well  known.  Several  extensions  to  the  theory  have 
occurred  more  recently,  made  possible  by  the  development  of 
modern  computers.  One  such  method  is  due  to  Mokry  [14]  who 
determined  the  corrections  for  subsonic  flow  in  tunnels  with 
perforated  walls.  Williams  and  Parkinson  [2],  in  studying 
the  effect  of  a  slotted  upper  wall,  used  a  singularity 
method  to  model  the  tunnel  flow.  De  Vries  and  Schipholt 
[15],  also  using  the  singularity  method,  found  large 
deviations  of  the  standard  forms  for  two  element  sections, 
and  showed  the  importance  of  the  model  rotation  point.  These 
modern  correction  techniques  are  not  in  wide-spread  use, 
probably  due  to  the  difficulty  in  setting  up  universal 
correction  formulae  for  non-standard  tunnel  geometries. 


F.  Conclusion 

The  main  result  of  this  brief  survey  of  wind  tunnel 
correction  theory  is  that  the  usual  corrections  are  valid 
for  small,  thin  airfoil  sections,  operating  at  low 
incidence,  at  or  near  the  tunnel  center-line.  When  all 
higher  order  terms  in  thickness,  chord-to-height  ratio,  and 
incidence  are  neglected,  the  lift  interference  is 
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independent  of  the  blockage,  and  vice  versa.  Clearly, 
however,  the  strict  assumptions  of  these  theories  are  not 
met  when  testing  large,  high-lift  models.  In  particular,  the 
treatment  of  multi-element,  high-lift  sections  with  the 
extreme  downward  deflection  of  the  trailing  streamlines  is 
not  complete.  Using  an  improved  method,  the  coupled  effects 
of  lift  interference  and  blockage  may  be  obtained 
simultaneously.  The  proposed  method,  described  in  the 
following  chapter,  utilizes  a  numerical  solution  for  the 
complete  flow  field  within  the  tunnel,  following  the 
suggestion  of  Garner,  without  restrictions  on  the  tunnel 
geometry . 


III.  Theoretical  Analysis  of  Tunnel  Interference 


This  chapter  will  describe  the  method  developed  to 
predict  the  combined  solid  blockage  and  streamline  curvature 
corrections  for  two  dimensional  testing  in  closed  wind 
tunnels.  These  corrections  are  necessary  due  to  the 
reduction  in  tunnel  cross  sectional  area  at  the  model  and 
the  constraining  effect  of  the  tunnel  walls,  and  may  be 
derived  using  potential  flow  analysis.  The  method  is  based 
on  the  theoretical  calculation  of  the  section  force  and 
moment  coefficients  for  both  a  free  field  and  a  bounded 
field.  The  difference  in  values  corresponds  to  the 
corrections  required  to  the  measured  wind  tunnel  data  to 
compensate  for  tunnel  wall  interference. 

Solutions  for  incompressible  potential  flow  in  an 
unbounded  fluid  have  generally  been  accomplished  using 
either  conformal  transformation  or  singularity  methods.  In 
the  first  method,  the  field  is  transformed  to  one  of  simple 
geometry  for  which  the  flow  solution  is  known,  usually 
uniform  flow  past  a  circular  cylinder  with  circulation.  A 
description  of  this  procedure  may,  be  found  in  the  usual 
major  references,  for  example  Abbott  and  von  Doenhoff  [16] 
or  Milne-Thomson  [17].  Examples  are  the  Joukowski, 
Karman-Tref f tz ,  and  von  Mises  airfoils  which  can  be  mapped 
directly  onto  the  circle.  These  mappings  are  restricted  to 
particular  families  of  airfoils,  few  of  which  are 
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particularly  useful  in  practice.  The  method  of  Theodorsen 
uses  a  double  mapping  for  an  arbitrary  airfoil,  first  to  a 
near-circle,  with  the  second  mapping  resulting  in  a  perfect 
circle.  A  conformal  transformation  technique  due  to  Williams 
[18]  allows  the  analysis  of  two  element  airfoil  sections  by 
mapping  to  two  circles  usinq  a  two  step  mapping  process. 

The  singularity,  or  panel,  methods  rely  on  the 
principle  of  superposition  of  solutions  for  the  linear 
Laplace  equation  for  the  velocity  potential  or  stream 
function.  The  airfoil  is  represented  by  a  distribution  of 
source  and/or  vortex  elements,  either  on  the  airfoil  surface 
or  in  the  interior  of  the  profile.  The  singularity  strengths 
are  found  such  that  the  total  flow  composed  of  the  onset 
uniform  stream  and  the  superposed  singularities  satisfies 
certain  boundary  conditions  on  the  airfoil  surface  and  at 
the  trailing  edge.  This  results  in  the  solution  of  a  large 
set  of  simultaneous  linear  algebraic  equations.  These 
integral  techniques  have  reached  a  high  level  of  refinement, 
can  be  used  for  multiple  bodies  of  arbitrary  shape,  and  can 
be  extended  to  three  dimensional  flows.  Both  de  Vries  and 
Schipholt  [15]  and  Williams  and  Parkinson  [2]  have  applied 
the  singularity  method  to  the  problem  of  wind  tunnel 
corrections.  Singularity  techniques  do  not  easily  provide 
the  flow  solution  throughout  the  field,  but  rather  give  the 
flow  on  the  airfoil  surface  as  primary  solution. 

The  method  to  be  used  in  the  present  work  is  quite 
different  from  these.  The  potential  flow  is  found  throughout 
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the  wind  tunnel  flow  field  using  finite  difference  methods 
to  solve  the  Laplace  equation  for  the  stream  function.  The 
influence  of  the  tunnel  geometry  is  felt  through  the 
boundary  conditions  imposed  on  the  solution  at  the  airfoil 
surface  and  on  the  tunnel  walls.  One  of  the  prime 
difficulties  with  finite  difference  techniques  is  the 
numerical  implementation  of  boundary  conditions  at 
irregularly  shaped  boundaries.  The  boundary  conditions  have 
a  dominant  effect  on  the  solution  since  the  Laplace  equation 
is  elliptic.  It  is  therefore  essential  that  the  boundary 
conditions  be  implemented  accurately  in  the  numerical 
scheme.  In  a  Cartesian  field,  special  finite  difference 
expressions  must  be  developed  for  nodes  on  or  near  an 
irregular  boundary,  using  special  differencing  and 
interpolation  techniques.  This  not  only  increases  the 
numerical  error  but  also  adds  considerably  to  the 
programming  effort  and  precludes  the  generality  of  the 
program  code  to  other  shapes. 

The  use  of  natural  coordinate  systems,  in  which  a 
boundary  contour  coincides  with  a  constant  coordinate  line, 
alleviates  this  difficulty  by  allowing  the  placement  of 
computational  nodes  exactly  on  the  flow  boundaries  and 
distributing  the  nodes  in  a  rational  manner  throughout  the 
field.  This  involves  a  transformation  of  coordinates,  much 
as  occurs  in  many  analytical  solutions  where  non-Cartesian 
coordinate  systems  are  chosen  to  facilitate  the  solution. 

The  governing  equations  are  written  in  the  transformed 
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coordinates,  and  boundary  conditions  are  carried  over  to  the 
transformed  plane.  The  solution  in  the  transformed  plane  may 
then  be  carried  out  on  a  regular  rectilinear  grid,  with  well 
posed  boundary  conditions.  A  simple  example  of  such  a 
coordinate  transformation  is  the  use  of  polar  coordinates, 
which  is  useful  when  the  physical  boundaries  are  composed  of 
concentric  circles,  as  shown  in  Figure  4a.  The  actual 
solution  of  the  transformed  governing  equation  takes  place 
in  the  transformed  plane.  Figure  4b.  Whether  the  solution  is 
analytic  or  numerical  in  nature,  the  use  of  the  coordinate 
transformation  has  allowed  the  boundary  conditions  to  be 
specified  on  straight  lines.  Of  particular  importance  to 
numerical  solutions  is  the  fact  that  the  flow  field  is  now 
represented  by  a  grid  of  equally  spaced  computational  nodes 
arranged  on  a  simple  rectangular  grid. 

There  are  many  possible  transformations  which  result  in 
such  rectangular  grids  and  in  general  such  curvilinear 
coordinate  systems  are  non-orthogonal .  These  coordinate 
systems  may  be  generated  in  several  ways,  which  may  be 
classified  as  algebraic,  conformal,  analytic,  or  numerical 
schemes.  The  application  of  analytic  and  conformal  grid 
generation  is  limited  in  that,  generally,  the  variety  of 
shapes  that  can  be  treated  is  small,  and  only  a  few 
specialized  mappings  can  handle  multiple  bodies.  However, 
analytic  and  conformal  mapping  techniques  generally  produce 
orthogonal  grids,  which  in  a  numerical  solution  reduces  some 
of  the  truncation  errors  which  result  in  the  discretized 
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Figure  4  Analytical  Coordinate  Transformation  -  Polar 
Coordinates 
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system.  Algebraic  systems  are  currently  under  development, 
see  for  example  Eiseman  [19]. 

The  method  of  coordinate  transformation  to  be  used  here 
is  a  numerical  grid  generation  technique.  This  method  is 
much  more  general  in  that  it  is  possible  to  generate  a 
coordinate  grid  for  a  flow  field  with  irregular  outer 
boundaries  containing  any  number  of  arbitrarily  shaped 
bodies.  The  grid  generation  and  flow  solution  algorithms  for 
different  flow  problems  require  only  that  the  flow  fields  be 
topologically  equivalent.  The  numerical  grid  generation 
methods  are  also  easily  extended  to  three  dimensional 
problems.  Although  this  method  of  grid  generation  was 
introduced  in  the  mid-1960’s,  its  recent  popularity  is  due 
largely  to  the  work  of  Thompson,  Thames,  Mastin,  et  al.,  at 
Mississippi  State  University  in  the  past  six  years. 
Application  of  the  method  of  boundary  fitted  coordinate 
systems  is  currently  being  made  to  a  wide  range  of  fluid 
flow  problems,  including  time-dependent,  compressible 
viscous  flow  with  movable  boundaries.  There  is  a  growing 
literature  on  the  use  of  these  coordinate  systems,  with 
recent  efforts  directed  towards  the  adaptive  and 
time-dependent  grid  systems  in  which  the  coordinate  grid  is 
adjusted  during  the  flow  solution  to  concentrate  the 
computational  nodes  in  regions  where  the  solution  exhibits 
large  gradients,  thus  increasing  the  solution  accuracy. 

In  the  method  of  Thompson  [20,21],  the  curvilinear 
coordinate  grid  is  generated  by  the  numerical  solution  of  a 
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set  of  elliptic  equations,  called  the  grid  generators. 
Boundary  conditions  for  the  grid  generators  are  specified  by 
the  desired  distribution  of  computational  nodes  along  the 
physical  flow  boundaries.  The  solution  of  the  generating 
system  establishes  a  correspondence,  or  mapping,  of  discrete 
points  in  the  physical  plane  to  the  nodes  in  the  rectangular 
computational  plane.  By  transforming  the  grid  generation 
equations  to  the  computational  plane,  the  generator  solution 
is  itself  carried  out  in  a  manner  which  takes  advantage  of 
the  well  posed  boundary  conditions. 

Once  the  coordinate  system  has  been  obtained,  the  task 
remains  to  solve  the  partial  differential  equation  which 
applies  to  the  particular  problem.  The  governing  equation  is 
transformed  to  the  new  coordinate  system,  just  as  is  the 
case  when  using  polar  coordinates  in  an  analytical  solution 
as  discussed  previously.  In  the  present  case,  the  potential 
flow  solution  is  found  by  transforming  the  Laplace  equation 
for  the  stream  function  and  solving  the  resulting  equation 
numerically  using  appropriate  boundary  conditions.  The 
stream  function  values  are  specified  on  the  flow  boundaries, 
i.e.  the  airfoil  surface  and  the  tunnel  walls,  entrance,  and 
exit  planes.  By  virtue  of  the  coordinate  transformation,  the 
boundary  computational  nodes  fall  on  rectilinear  boundaries 
in  the  computational  grid. 

For  lifting  potental  flow,  the  total  solution  is  found 
by  superposition  of  the  basic  solutions  for  uniform  onset 
flow  with  zero  circulation  and  for  zero  onset  flow  with  unit 
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circulation,  which  are  obtained  separately.  The  required 
circulation  is  determined  by  application  of  one  of  several 
forms  of  the  Kutta  condition  which  are  also  expressed  in  the 
computational  plane.  Since  each  computational  node  in  the 
transformed  plane  corresponds  uniquely  to  some  point  in  the 
physical  plane,  the  total  flow  solution  is  therefore  known 
at  these  points. 

The  remainder  of  the  analysis  is  composed  of  the 
calculation  of  the  surface  pressures  on  the  airfoil  and 
tunnel  walls,  which  may  be  integrated  to  give  the  lift  and 
pitching  moment  coefficients.  All  of  these  calculations  may 
be  carried  out  in  the  computational  plane,  and  in  fact,  once 
the  physical  boundary  points  are  input  into  the  grid 
generation,  they  never  re-appear  until  the  final  results  are 
presented. 

The  preceeding  method  of  analysis  was  used  to  obtain 
the  flow  solution  for  an  airfoil  section  in  a 
two-dimensional  wind  tunnel.  Free  field  values  for  the  force 
and  moment  coefficients  were  obtained  using  the  singularity 
method  of  Kennedy  [13].  By  comparing  results  from  these  two 
analyses  the  required  adjustments  needed  to  correct  the  wind 
tunnel  data  for  the  interference  of  the  walls  was  obtained. 

A.  Grid  Generation 

In  the  method  of  Thompson,  the  transformation  to  ( ^  ,  yj  ) 
coordinates  represented  by  ^(x,y)  and  ^(x, y),  is  found  by 
the  solution  of  two  elliptic  partial  differential  equations 
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(the  grid  generators): 


(22a) 


72,1  =  nxx  +  nyy  *  ^.n) 


(22b) 


with  Dirichlet  boundary  conditions.  The  boundary  conditions 
take  the  form  of  the  desired  coordinate  values  (£,n)  at 
discrete  points  on  the  irregular  physical  boundaries.  One 
coordinate  is  held  constant  while  the  other  varies 
monotonically  around  or  along  the  boundary.  This  form  of 
grid  generator  and  boundary  conditions  ensures  a  unique 
mapping  from  the  physical  to  the  computational  space.  The 
inhomogeneous  functions  P(£,n)  and  Q(£,n)allow  the  generated 
mesh  to  be  concentrated  in  regions  of  interest,  usually 
where  gradients  in  the  flow  solution  are  high  or  the 
boundaries  are  highly  curved. 

Since  it  is  desired  to  perform  all  calculations  in  the 
regular  transformed  plane,  i.e.  to  find  x(£,n)and  y(£#n),  the 
dependent  and  independent  variables  in  the  above  generating 
equations  (21)  are  interchanged  to  give: 


ax  -  2Sx  +  yx  =  -  J2 

££  nn 


Px_  +  Qx 

K  n 


(23a) 


(23b) 


where 
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2  2 

a  =  x  +  y 
n  n 


3  -  x5xn  +  y?yn 


Y 


with  boundary  conditions  (i.e.  values  of  x  and  y  at  the 
desired  computation  points  oh  the  physical  boundaries) 
specified  on  the  straight  boundaries  of  the  computational 
grid. 

The  terms  a  ,  3  ,  y  ,  and  J  (the  Jacobian  of  the 

transformation)  are  called  the  transformation  metrics  or 
scale  factors.  The  metrics  are  obtained  during  the  grid 
generation,  i.e.  the  solution  of  equations  (23a)  and  (23b), 
and  will  appear  later  in  the  flow  solution  equations.  A 
physical  interpretation  of  the  metrics  may  be  made  in  terms 
of  the  grid  orientation,  orthogonality,  cell  size,  and  cell 
aspect  ratio,  as  discussed  in  a  later  section. 

There  are  numerous  configurations  which  may  be 
generated  using  the  elliptic  system  due  to  the  flexibility 
allowed  in  the  specification  of  the  boundary  conditions  on 
the  computational  grid.  For  the  present  work,  four  grid 
types  were  investigated.  The  first  two  are  commonly  referred 
to  in  the  literature  as  O-type  and  C-type  grids.  The  last 
two  are  designated  here  as  S-type  and  2S~type  grids. 


Grid  Types 

The  O-type  grid  as  it  appears  in  the  physical  plane  is 
shown  in  Figure  5a.  The  corresponding  computational  grid  is 
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Figure  5  0-Type  Grid  -  General  Form 
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shown  in  Figure  5b,  with  numbered  points  corresponding  to 
those  in  the  figure  above.  The  notable  features  of  this  grid 
type  are  the  uniformity  of  the  mesh  near  the  inner  boundary, 
and  the  re-entrant  line  on  the  left  and  right-hand  sides  of 
the  computational  grid  which  appears  as  a  cut  in  the 
physical  plane  between  the  inner  and  outer  boundaries.  The 
O-type  .grid  is  similar  to  a  polar  coordinate  system,  where 
the  re-entrant  line  occurs  at  9  =  0  and  2tt  .  Solutions  are 
required  to  be  continuous  across  the  re-entrant  line,  and  in 
fact  the  solution  is  periodic  in  .  Values  of  x  and  y  are 
neither  specified  nor  allowed  along  the  re-entrant  line,  but 
are  found  during  the  grid  generation  by  requiring  continuity 
of  the  solution  and  first  derivatives  across  this  line,  as 
shown  by  the  symbols. 

The  C-type  grid  is  shown  in  Figure  6a  and  6b,  again 
with  numbered  points  to  establish  the  correspondence  of 
computational-  points  with  the  physical  plane.  This  grid  also 
has  a  re-entrant  line,  this  time  along  a  portion  of  the 
lower  edge  of  the  computational  grid.  Again,  specification 
of  values  of  x  and  y  at  these  points  is  neither  necessary 
nor  allowed. 

The  S-type  and  2S-type  grids  present  another  possible 
arrangement  in  the  transformed  plane.  For  these  grids, 
interior  boundaries  are  mapped  onto  slits  in  the  transformed 
plane.  For  a  single  interior  body,  this  mapping  results  in  a 
grid  as  shown  in  Figure  7a,  with  the  slit  shown  in  the 
transformed  plane,  Figure  7b.  For  the  grid  generation,  the 
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Figure  6  C-Type  Grid  -  General  Form 
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b)  transformed  plane 


Figure  7  S-Type  Grid  -  General  Form 
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slit  boundary  condition  is  multi-valued  in  x  and  y,  the 
value  used  depending  on  which  side  of  the  slit  computations 
are  performed.  There  is  no  re-entrant  line.  The  25-type,  or 
double  slit  geometry,  is  shown  in  Figure  8a  and  8b  for  two 
interior  bodies.  Here  also  the  slit  boundary  conditions  are 
multi-valued.  Using  the  slit  geometries,  it  is  possible  to 
set  up  grids  for  much  more  complicated  sections:  a  4S-type 
grid  might  be  used  for  an  airfoil  section  with  a  double 
slotted  flap  and  leading  edge  slat.  The  slits  may  also 
appear  vertically  in  the  transformed  plane.  The  present  work 
has  been  limited  to  the  simplest  two  element  case  of  a 
single  slotted  flap. 

Of  these  four  grid  types,  only  the  O-type  grid  has  been 
investigated  in  detail  for  potential  flow  solutions 
(Thompson [ 24 ]) .  Each  of  these  grid  configurations  offered 
some  advantages  to  some  part  of  the  problem  at  hand,  and 
were  therefore  investigated  in  some  detail.  The  over-all 
accuracy  of  the  flow  solution  depends  on  several  properties 
of  the  generated  coordinate  grid  in  the  physical  plane.  In 
general,  as  with  all  finite  difference  methods,  accuracy 
deteriorates  if  the  grid  cells  are  large,  or  if  the  grid 
geometry  changes  abruptly  or  rapidly  from  point  to  point  in 
the  physical  plane.  This  is  especially  true  if  the  solution 
exhibits  high  gradients  in  these  regions  of  rapid  grid 
change.  Furthermore,  although  it  is  not  necessary  that  the 
grid  be  orthogonal,  highly  skewed  grids  contribute  to  the 
local  truncation  error  of  the  flow  solution.  The  effect  of 
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Figure  8  2S-Type  Grid  -  General  Form 
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the  physical  geometry  of  the  computational  grid  is  felt 
through  the  metrics.  Since  the  metrics  are  calculated  using 
standard  finite  difference  forms,  they  become  susceptible  to 
the  normal  errors  of  truncation.  Thus,  although  the 
transformed  flow  equation  may  be  exact  in  form,  the  finite 
difference  solution  will  be  degraded  by  truncation  errors 
arising  not  only  from  the  discretization,  but  also  from  the 
transformation  metrics  which  it  utilizes,  as  demonstrated  by 
Kerlick  and  Klopfer  [22].  It  is  desireable,  therefore,  to 
obtain  grids  with  near  uniform  cell  size,  or  no  more  than 
slow  size  variation  over  the  grid,  and  with  no  more  than 
slight  skewness.  With  this  in  mind,  representative  grids  of 
each  type  for  the  wind  tunnel  geometry  are  shown  in  Figures 
9  to  13. 

The  O-type  grid.  Figure  9,  has  the  best  mesh  properties 
near  the  airfoil  surface,  with  the  exception  of  the  sharp 
trailing  edge,  but  very  poor  properties  near  the  outer 
boundaries,  where  the  mesh  is  composed  of  large,  highly 
skewed  cells.  In  these  outer  regions  of  the  flow  the 
solution  does  not  change  rapidly,  so  the  grid  may  be 
tolerated. 

The  C-type  grids,  Figure  10,  were  investigated  as  a 
means  for  improving  the  grid  quality  upstream  from  the 
model.  This  was  desired  in  order  that  the  wind  tunnel 
contraction  could  be  included  in  the  solution  field  if 
desired,  as  shown  in  Figure  11. 
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b)  airfoil  detail 


gure  9  0-Type  Grid  -  Wind  Tunnel 
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a)  entire  flow  field 


b)  airfoil  detail 


Figure  10  C-Type  Grid  -  Wind  Tunnel  Geometry 


■ 


54 


Figure  1 1  C-Type  Grid  -  Tunnel  Contraction  Geometry 
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The  S-type  grids,  which  are  well  suited  to  duct  flows 
as  they  have  excellent  form  upstream  and  downstream  from  the 
model,  are  poor  near  the  section  leading  edge  as  shown  in 
Figures  12  and  13.  The  S-type  grids  are  also  well  suited  to 
the  contraction  geometry.  The  two  element  section  grids 
could  be  generated  using  a  form  of  the  C-type  grid,  however 
the  mesh  generated  in  the  slot  region  is  unsatisfactory  due 
to  extreme  skewness  and  poor  spacing.  Thompson  solved  this 
problem  (see  [21])  by  performing  an  auxilliary 
transformation  where  the  section  coordinates  were  inverted 
with  respect  to  infinity.  This  had  the  added  advantage  of 
allowing  the  outer  boundary  to  be  moved  a  very  large 
distance  away  from  the  airfoil,  and  the  implementation  of 
near-infinity  boundary  conditions  for  free  field  flow 
solutions.  The  added  complexity  in  the  flow  equations,  which 
must  also  be  inverted,  prevented  this  method  from  being 
applied  for  an  actual  flow  solution  [23].  Using  the  2S-type 
grids  an  adequate  mesh  could  be  generated  in  the  slot  region 
with  little  difficulty. 

It  is  observed  that  none  of  these  grids  is  uniformly 
orthogonal,  however  this  property  is  not  required,  nor  even 
desired.  In  fact,  demanding  an  orthogonal  grid  would 
over-specify  the  problem  and  eliminate  the  freedom  with 
which  the  boundary  point  distribution  could  be  specified. 

Minimization  of  the  grid  induced  errors  may  be 
accomplished  either  by  proper  specification  of  the  boundary 
point  distribution,  or  by  use  of  the  weighting  functions  P 
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a)  entire  flow  field 


b)  airfoil  detail 


Figure  12  S-Type  Grid  -  Wind  Tunnel  Geometry 
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a)  entire  flow  field 


b)  airfoil  detail 


Figure  13  2S-Type  Grid  -  Wind  Tunnel  Geometry 
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and  Q. 


Boundary  Point  Distribution 

The  over-all  accuracy  is  dependent  on  the  chosen 
distribution  of  points  on  the  boundary  contours.  It  is 
desirable  to  use  as  few  points  as  possible,  consistent  with 
desired  accuracy,  in  order  to  reduce  required  computer  time. 
Critical  regions  of  the  flow  are  near  the  section  leading 
edge  where  the  flow  curvature  is  high,  and  the  trailing  edge 
where  the  Rutta  condition  is  implemented.  Thompson  [24] 
showed  that  for  the  O-type  grids  accuracy  was  improved  more 
by  concentration  of  points  near  the  trailing  edge  than  near 
the  leading  edge.  In  the  present  work  this  was  also  found  to 
be  the  case  for  the  other  grid  types. 

The  airfoil  boundary  points  were  specified  using  the 
distribution : 

Xi  =  \  (1  •  COS(^i) 

where  -ji 

♦i  =  nFT  %  ’  1  =  1  ’  NF 

4>i  *  %  +  I^xrr  »  1  =  NF  +  1  ,  NF  +  NA 

giving  NF  points  distributed  between  the  leading  edge  (with 
slight  concentration)  and  the  point  corresponding  to  4>^p  , 
usually  120  degrees,  and  NA  points  between  <f>^p  and  the 
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trailing  edge.  With  these  x  coordinates,  the  corresponding  y 
coordinates  of  the  surface  were  found  using  a  cubic  spline 
interpolation  routine,  following  the  procedure  of  Kennedy 
[13],  using  the  section  input  coordinates.  A  typical  point 
distribution  is  shown  in  Figure  14  for  a  20%  thick 
symmetrical  Joukowski  airfoil,  using  NF=NA=15. 

The  outer  bpundary  point  distribution  scheme  was 
identical  for  each  of  the  four  grid  types  investigated.  The 
cross  stream  boundary  points  were  equally  spaced  along  the 
tunnel  inlet  and  exit  planes,  while  the  wall  point  spacing 
increased  from  the  region  of  the  model  to  the  end  lines. 

Grid  Control  Functions 

In  order  to  optimize  the  grid  within  the  field,  a  high 
degree  of  control  is  possible  using  the  inhomogeneous 
functions  P  and  Q  in  the  grid  generators.  By  setting  the 
magnitude  and  sign  of  the  functions  P  and  Q  at  each  point, 
the  generated  grid  can  be  tailored  to  the  problem  at  hand. 
The  form  of  the  control  functions  used  in  the  present  work 
is  the  same  as  that  of  Thompson  [21].  The  functions  are 
evaluated  at  each  point  in  the  transformed  plane  as  sums  of 
decaying  exponentials,  where  the  magnitude  and  sign  of  each 
term  are  arranged  to  attract  or  repel  grid  lines  from 
specified  locations.  The  strength  of  the  attraction  and 
eventual  clustering  of  cells  is  determined  by  the  magnitude 
of  the  function,  while  the  sign  determines  the  direction  of 
attraction  of  the  coordinate  lines  in  the  physical  plane. 
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+  boundary  points  used  for 
grid  generation 


Figure  14  Airfoil  Boundary  Point  Distribution 
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The 


functional 

Q(£>n)  = 


forms  of  P  and  Q  are: 
n  r 

-  I  a.  sgn(n-nJ  exp  - 
1*1  L 


m 

-  y  b.  sgn(n-n.)  exp 
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(24a) 


P(5»n) 
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sgn(^-C1- )  exp 


where 


m 

•  s  b.  sgn(5-£.)  exp 
j=l  J  J 


(24b) 
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The  number  of  terms  in  each  summation  -is  determined  by  the 
desired  complexity  of  the  grid  control  desired.  The  ,bj 
and  c-,d-  are  positive  amplitude  and  decay  factors, 

"  v 

respectively,  which  are  not  necessarily  the  same  for  both 
functions.  The  presence  of  the  "sgn"  function  in  (24a)  and 
(24b)  allows  the  sign  of  the  control  function  to  change  on 
either  side  of  the  desired  location  of  concentration. 
Positive  values  for  P  and  Q  tend  to  move  grid  lines  towards 
higher  values  of  £,  and  ,  respectively,  and  vice  versa. 
Thus,  the  first  summation  in  equation  (24a)  results  in 
attraction  of  77  -coordinate  lines  toward  the  given  77-  line 
in  the  physical  plane  while  the  second  summation  causes 
attraction  of  77  lines  to  a  specific  point  •  Similar 
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results  are  obtained  with  equation  (24b),  where  now  the  E, 
lines  are  affected.  Examples  of  grid  control  using  the 
weighting  functions  are  shown  in  Figures  15  to  18,  using  the 
same  boundary  point  distributions  as  Figures  9  to  13. 

Overall  grid  size  and  control  amplitudes  and  decay  factors 
used  are  indicated  below  each  figure.  In  Figure  15,  rj  lines 
have  been  concentrated  near  both  the  outer  boundary  and  the 
airfoil  surface  by  using  the  line  attraction  form  of  the  Q 
control  function,  the  first  term  in  (24a),  to  attract 
lines  to  both  the  upper  and  the  lower  boundary  of  an  O-type 
grid.  In  Figure  16,  the  P  function  has  been  used  to  attract 
£  lines  toward  the  leading  and  trailing  edge  points,  the 
second  term  in  (24b),  using  an  S-type  grid.  For  the  S-type 
grids  in  the  wind  tunnel  geometry,  it  was  found  to  be 
necessary  to  use  a  two-step  decay  function  to  prevent  the 
E,  lines  near  the  tunnel  end  lines  from  responding  to  the 
high  control  values  required  near  the  leading  and  trailing 
edges,  as  shown  upstream  of  the  model  in  Figure  16b.  An 
example  of  a  2S-type  grid  generated  with  and  without  extreme 
grid  control  is  shown  in  Figures  17  and  18.  In  Figure  18, 
both  line  and  point  attractions  for  P  and  Q  were  necessary, 
using  several  terms  in  the  summations  of  equations  (24a)  and 
(24b),  to  obtain  an  adequate  grid  distribution  near  the 
model . 

Guidelines  for  the  use  of  grid  control  may  be  deduced 
from  the  tendency  for  the  Laplacian  operator  used  to 
generate  the  grids  to  distribute  the  grid  points  evenly  in 
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a)  no  grid  control 


b)  with  grid  control 

Tj  - 1 1  ns  attraction  to  inner  and  outer  boundaries 
amp  1 i tude  =  200 . 0  decay  rate=0.5 


Figure  15  Grid  Control  -  0-Type  Grid  with  Q  Control 
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a)  no  grid  control 


b)  with  grid  control 


point  attraction  of  £, -lines  to  airfoil  LE  and  TE 
amp  1 i tude  =  300 . 0  decay  rate=0.35 


Figure  16  Grid  Control  -  S-Type  Grid  with  P  Control 
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Figure  17  2S-Type  Grid  -  Airfoil  Detail  -  No  Grid  Control 
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Grid  Control 
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the  field  and  to  avoid  concave  boundaries  and  move  toward 
convex  boundaries.  In  order  to  maintain  the  desired  grid 
properties  extensive  use  of  the  grid  control  functions  is 
mandatory  to  achieve  accurate  results  for  the  flow  solution, 
especially  for  the  slit  geometries.  Some  very  interesting 
work  is  being  done  in  this  area,  see  for  example  Brackbill 
[25]  or  Anderson  and  Rai  [26],  where  the  grid  control  is 
implemented  to  minimize  some  form  of  a  penalty  function,  for 
example  the  non-orthogonality  of  the  grid,  or  to  respond  to 
high  gradients  in  the  flow  solution.  These  are  termed 
adaptive  grid  methods.  For  the  present,  the  optimal  grid 
control  parameters  were  determined  by  trial  and  error  and 
were  fixed  throughout  the  solution.  These  parameters  depend 
primarily  on  the  grid  type  and  size.  For  the'  grids  used  in 
the  present  work  (normally  100  X  80),  amplitudes  ranged  from 
100  to  600,  and  decay  factors  ranged  from  0.1  to  0.5. 

B.  Potential  Flow  Solution 

The  stream  function  for  the  two  dimensional,  inviscid, 
incompressible  flow  satisfies  Laplace's  equation 

V^iii  =  ib  +  ip  =  0 

v  ax  yyy 

with  Y  having  constant  values  on  solid  boundaries  (i.e. 
these  boundaries  must  be  streamlines),  and 
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ip  =  ycos0  -  xsinQ 

at  a  large  distance  from  a  body  immersed  in  a  uniform  flow 
at  angle  9  to  the  x-axis.  When  transformed  to  the 
computational  plane,  equation  (25)  becomes 

-  26^n  +  Y^nn  =  -  J2[P^  +  Q^]  (26) 

where  ot  ,  (3  ,  T  ,  and  J  are  the  transformation  metrics 

obtained  during  the  grid  generation  and  P  and  Q  are  the  grid 
control  functions.  The  boundary  conditions  in  the 
computational  grid  are  specified  along  lines  corresponding 
to  the  body  contours  and  the  outer  boundary  of  the  flow 
field.  In  the  present  case,  these  boundaries  are  either 
streamlines  of  the  flow  (  Y  =constant),  or  occur 
perpendicular  to  the  flow  direction  (  Y  a  linear  function 
of  position ) . 

Considering  first  a  single  element  section,  since  the 

Laplace  equation  is  linear,  the  total  flow  may  be  separated 

into  three  parts:  ^  due  to  a  uniform  onset  flow  at  zero 
2 

degrees,  ip  due  to  uniform  onset  . flow  at  ninety  degrees,  and 
3 

.ip  for  circulatory  flow  with  zero  onset  velocity.  These 
basic  solutions  are  found  individually  by  solving  equation 
(26)  with  the  boundary  conditions  shown  in  Table  1. 

The  total  flow  is  then  found  by  superposition  of  the 
three  base  solutions: 


69 


Table  1  Boundary  Conditions  for  Single  Element  Potential 
Flow 


Basic  Solution 

Unbounded 

ai  r foi 1 

F  low 

outer 

boundary 

Wind  Tunnel  Flow 

ai r  foi 1  outer 

boundary 

1 .  Uni  form  stream 

at  Zero  Degrees 

0 

y 

0  y 

2.  Uniform  stream 

at  Ninety  Degrees 

0 

-x 

“  •  “•  —  »  •*» 

3.  Circulatory  Flow 

1 

0 

1  0 

(27) 


T  1  ?  3 

ip  =  ip  cos6  +  ip  sinQ  +  XiJj 

where  9  is  the  desired  angle  of  the  onset  flow,  and  X  is 

to  be  determined  for  a  given  9  by  the  Kutta  condition.  For 

the  general  problem,  the  three  basic  solutions  need  to  be 

calculated  only  once.  Thereafter,  the  total  stream  function 

at  any  angle  may  be  found  by  finding  the  appropriate  X  ,  as 

discussed  in  the  following  section. 

For  the  two  element  case,  the  total  flow  is  separated 

12  3 

into  four  solutions:  ip  and  ip  as  before,  ip  for  circulatory 

4 

flow  about  the  main  section,  and  ip  for  circulatory  flow 
about  the  flap.  Again,  these  solutions  of  the  Laplace 
equation  are  found  separately  using  the  boundary  conditions 
given  in  Table  2.  The  total  flow  is  then  given  by 

^  =  ijJ  cosQ  +  ip^  sin0  +  (28) 

where  0  is  the  angle  of  the  onset  flow,  and  X^  and  X^ 
are  determined  by  the  satisfaction  of  the  Kutta  condition  on 
the  two  components  of  the  section.  This  involves  the 
simultaneous  solution  of  two  linear  algebraic  equations  for 
each  angle  of  attack,  as  discussed  in  the  following  section. 

For  the  present  problem,  the  onset  flow  and  tunnel 
walls  are  parallel,  and  the  model  incidence  is  set  by 
rotating  the  model  coordinates  about  a  given  point  on  the 
chord  line.  Hence,  9=0  ,  and  the  y1  term  vanishes  from 


. 


. 
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Table  2  Boundary  Conditions  for  Two  Element  Potential  Flow 


Unbounded  Flow  Wind  Tunnel  Flow 


Basic  Solution 


main  flap  outer  main  flap  outer 

boundary  boundary 


1 .  Uni  form  Stream 
at  Zero  Degrees 


2.  Uniform  Stream 

at  Ninety  Degrees  0  0  -x 


3.  Circulatory  Flow 

(main)  1 


0 


0 


4.  Circulatory  Flow 
(flap) 


-1  0 


1  0 
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the  total  flow.  For  the  single  element  section  this  leaves 
just  two  basic  flow  solutions  to  be  found  to  give  the  total 
flow  by 


(29) 


while  three  basic  solutions  are  required  for  two  element 
sections,  where  the  total  flow  is  given  by 


T  1  3  4 

ip  =  ip  +  A^ip  + 


(30) 


Boundary  conditions  for  these  basic  solutions  in  the  wind 
tunnel  geometry  are  also  given  in  Tables  1  and  2.  For  each 
new  geometric  angle  of  attack  of  the  section  in  the  wind 
tunnel,  a  new  grid  must  be  generated,  and  new  basic 
solutions  found,  which  involves  a  great  deal  of  extra 
computation.  If  the  change  in  angle  of  attack  is  not  large, 
the  basic  solutions  are  very  similar,  so  that  the  solution 
for  one  angle  may  be  used  as  input  to  the  iterative  solution 
at  the  new  angle. 


Evaluation  of  Velocities  and  Surface  Pressures 

Velocity  components  of  the  total  flow  solution  may  be 
found  by  evaluating 


u  =  ip 


T 


T 


y 


X 


9 


- 


. 
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which  in  the  computational  grid  take  the  form  (omitting  the 


superscript  "T"  for  clarity) 

u  =  !  [Vn  -  V?] 

(31a) 

v  -  j  [y?*n  -  yn*5] 

(31b) 

In  the  field,  these  derivatives  are  found  using  second  order 
central  difference  formulae.  On  the  boundaries,  second  order 
one-sided  expressions  are  used  for  the  derivatives  normal  to 
the  boundary.  Velocity  components  on  solid  boundaries  ( V 
=constant),  which  here  correspond  to  lines  of  constant 

are  given  by  (since  ^  =0) 

xrip 

__  u  =  ^  31  (32a) 


On  the  surface, 
by 


yrip 
=  £  n 


(32b) 


the  tangential  velocity  is  therefore  given 


2  o  y  ^ 

V  =  (u  +v  )  =  - - r-3 


The  pressure  coefficient,  either  in  the  field  or  on  the 
boundary  is  given  by 


Cp  ■  1  -  v2 

In  the  field,  then,  the  pressure  coefficient  is  given  by 


. 


.  ■ 


. 
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C  =  1  -  \ 
p  .l2 


a(*5)2  -  2Bip^n  +  Y(<Pn)2 


while  on  the  solid  boundaries 


y(\Y 

c  =  i  — 

P  ,2 


(33a) 


(33b) 


Kutta  Condition  Implementation 

The  amount  of  circulation  generated  by  a  lifting 
airfoil  in  an  inviscid  fluid  is  determined  by  the  condition 
that  the  flow  leave  the  trailing  edge  smoothly,  commonly 
referred  to  as  the  Kutta  condition.  Historically,  this 
condition  has  been  implemented  in  several  ways.  In  conformal 
transformation  methods,  the  circulation  is  given  around  the 
cylinder  in  the  transformed  plane  by  requiring  that  the  rear 
stagnation  point  map  onto  the  airfoil  trailing  edge.  In  the 
singularity  techniques,  the  Kutta  condition  has  been 
implemented  in  terms  of  velocities  or  flow  direction.  The 
velocity  condition  takes  two  forms:  equal  velocity  on  the 
upper  and  lower  surface  at  the  trailing  edge,  or  zero 
velocity.  A  less  restrictive  (and  physically  more 
reasonable)  implementation  is  specification  of  the  .flow 
direction  at  the  trailing  edge.  This  is  accomplished  (see 
Kennedy  [13])  by  placing  a  control  point  on  the  trailing 
edge  bisector  at  a  very  short  distance  behind  the  trailing 
edge  and  requiring  the  stream  function  at  the  control  point 
to  coincide  with  the  value  on  the  airfoil  surface.  An 


- 
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alternative  implementation  of  the  trailing  point  condition 
is  the  requirement  that  the  velocity  vector  at  the  trailing 
edge  be  in  the  direction  of  the  trailing  edge  bisector.  For 
the  singularity  techniques,  Kennedy  showed  the  trailing 
point  Kutta  condition  to  be  superior  to  the  velocity  methods 
by  allowing  a  reduction  in  the  number  of  airfoil  surface 
elements. 

Both  the  velocity  and  trailing  streamline  conditions 
were  investigated  for  each  grid  type  in  the  present 
numerical  scheme. 

The  equal  velocity  condition  is  implemented  as  follows 
for  a  single  element  section.  The  velocity  at  a  point  on  the 
surface  of  a  single  element  section  is  given  by 


V  = 


1/2.  T 

Y  ^ 

n 


1  2  3 

W  cos0  +  \p  sin©  +  Aiii 
n  n  n 


(34) 


where  the  are  found  using  one  sided  second  order  finite 
difference  formula.  Expressions  for  V,  containing  the 
parameter  A  are  found  for  points  near  the  trailing  edge  on 
the  upper  and  lower  surface.  Using  one,  two,  or  three 
neighbouring  points,  an  extrapolation  is  made  on  each 
surface  to  obtain  an  expression  for  the  velocity  at  the 
trailing  edge.  As  the  extrapolated  trailing  edge  velocity  is 
to  be  the  same  for  both  the  upper  and  lower  surface,  an 
equation  results  which  gives  A  explicitly  in  terms  of  the 
base  solution  derivatives  and  0  .  The  extrapolation  formula 


. 
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takes  one  of  the  forms  (Thompson  [24]) 


Vth-1  =  VTE  =  VTE+1 


(35a) 


2VTE-1  _  vTE-2  =  VTE  =  2VTE+1  .  vTE+2 


(35b) 


(35c) 


where  the  superscripts  indicate  the  position  relative  to  the 
trailing  edge  as  shown  in  Figure  19.  The  first  condition, 
equation  (35a),  results  in  a  stagnation  point  at  the 
trailing  edge,  the  zero  velocity  condition.  The  second 
condition,  equation  (35b),  represents  a  linear  extrapolation 
to  the  trailing  edge  velocity,  while  the  third  expression, 
(35c),  results  from  a  quadratic  extrapolation.  Thompson 
found  a  slight  improvement  in  accuracy  for  O-type  grids 
using  the  quadratic  scheme.  A  similar  result  occurs  for  the 
C-type  grids  due  to  the  identical  nature  of  the  grid  near 
the  trailing  edge  for  these  two  mappings.  For  the  S-type 
grids,  it  was  necessary  to  re-examine  this  problem,  where  a 
similar  improvement  was  found. 

For  the  2S-type  grids  the  parameters  X,  and  XL  are 
found  by  extrapolations  to  the  two  trailing  edge  velocities 
on  the  two  components.  This  results  in  equations  of  the  form 


' 
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b)  transformed  plane:  O-type,  C-type  grids 
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c)  transformed  plane:  S-type  grids 


Figure  19  Trailing  Edge  Points  For  Velocity  Kutta  Condition 
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(36a) 


(36b) 


M  F 

where  the  and  Cj  refer  to  main  and  flap  sections,  and  are 

‘  i 

formed  from  the  derivatives  of  the. basic  solutions  Y  ,  y  , 
and  Y  at  the  extrapolation  points  and  the  desired  flow 


angle  9  . 

The  trailing  streamline  condition  was  implemented  only 
for  the  0-  and  C-type  grids.  The  velocity  vector  at  a 
computational  node  (eg.  the  trailing  edge)  is  given  by 

V  =  uT  +  vj 

and  the  bisector  at  the  trailing  edge  is  given  by 

b  =  b-,1  +  b2J  ■ 

If  the  velocity  vector  and  bisector  coincide,  then 


V  x  b  =  0 


hence 


ub^  -  vb^  =  0 


(37) 


■ 


. 
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the  velocity  components  are  given  by  equation  (32).  Using 

y 

the  full  expression,  equation  (27),  for  ly  ,  equation  (37) 
reduces  to 

1  2 
cos@  +  ip  sine 

X  =  -  -n - J* -  (38) 

This  is  a  much  simpler  expression  than  the  equal  velocities 
method,  as  only  three  finite  difference  derivatives  need  to 
be  evaluated  at  the  trailing  edge.  Moreover,  the  grid  is 
most  dense  at  the  trailing  edge,  so  that  accuracy  is 
maintained.  For  the  S-type  grids  a  similar  expression  does 
not  result,  due  to  the  multi-valued  coordinate  mapping  on 
the  slit  at  the  trailing  edge,  and  recourse  must  be  made  to 
the  velocity  condition. 

The  implementation  of  the  Kutta  condition  was  found  to 
be  the  most  critical  aspect  of  the  present  work.  None  of  the 
methods  described  above  were  found  to  be  applicable  on  all 
four  grid  types.  In  general,  the  velocity  conditions  worked 
best  in  the  S-  and  2S-type  grids,  while  the  trailing  point 
condition  worked  best  for  the  0-  and  C-type  grids. 


Lift  and  Moment  Calculations 

The  force  system  on  the  airfoil  is  shown  in  Figure  20. 
The  chordwise  force  coefficient  Cx  is  negative  for  positive 
angles  of  attack  in  an  unbounded,  inviscid  flow  so  that  the 
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Figure  20  Airfoil  Force  System 


. 
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resultant  force,  Cr  ,  is  perpendicular  to  the  free  stream. 
The  force  and  pitching  moment  coefficients  are  found  by 
integration  of  the  pressure  distribution,  which  can  be  done 
in  several  ways.  In  the  present  work,  these  integrations 
were  performed  in  the  physical  plane  using  the  relations 
(see  Houghton  and  Brock  [27]) 


CN  = 


CPS.  '  CPu 


dx 


y 

/max  r 


V 


y  . 
min 


Cp  -  Cp 

■  f  -* 


dy 


,1  r— 


'M 


LE 


CPe  ’  CPu 


x  dx 


(39a) 

(39b) 


(39c) 


to  obtain  normal  and  chordwise  force  coefficients,  and  the 
leading  edge  pitching  moment  coefficient.  The  contribution 
to  the  leading  edge  pitching  moment  due  to  chordwise  forces 
was  ignored,  as  it  is  very  small.  The  integrations  were 
carried  out  using  the  trapezoidal  rule.  The  force  components 
were  resolved  to  give  the  lift,  drag,  and  quarter  chord 
pitching  moment  coefficients  using 


C£  =  cose  -  C  sine 


(40a) 


■ 
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C  ,  =  CM  si n0  +  C  v  cos0 


(40b) 


(40c) 


These  coefficients  may  also  be  evaluated  directly  in 
the  transformed  plane  using  the  expressions 

N 


(41a) 


C  =  -  C„  xrcos0  +  yrsin0  d£ 
&  P  £  £ 


N 


(41b) 


C 


d 


where  the  integrals  may  be  easily  evaluated  using  either  the 
trapezoidal  or  Simpson's  rule  over  the  N  coordinate  points 
which  define  the  airfoil  surface.  Thompson  [24]  used  this 
procedure  for  O-type  grids,  noting  that  accuracy 
deteriorated  using  Simpson's  rule. 

The  resultant  force  coefficient  may  be  obtained  by  the 
calculation  of  the  circulation,  F  ,  around  the  section 
which  is  defined  by 


r  =  -  o  v  •  ds 


where  the  circulation  has  been  normalized  with  respect  to 


■ 
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Vo  .  For  the  O-type  grid  the  line  integral  is  evaluated  in 
the  transformed  plane  on  a  line  of  constant  rj  which  is  the 
map  of  a  closed  curve  in  the  physical  plane.  The  integral 
becomes 

N 

►  f 

(I  V  •  d?  =  - 
*  • 

All  derivatives  in  the  integrand  are  evaluated  by  second 
order  central  difference  expressions  on  the  7}  -line  nearest 
the  airfoil  surface  and  the  integral  is  evaluated  with  the 
trapezoidal  rule  using  Romberg  improvement.  Other 
integration  schemes  (Simpson’s  rule  and  spline  integration) 
were  investigated,  with  no  appreciable  improvement  in 
accuracy.  The  force  coefficient,  Cr  ,  is  found  using 


This  procedure  provides  another  check  on  the  lift 
calculation,  but  was  not  implemented  in  the  C-  or  S-type 
grids . 

Both  the  force  and  pitching  moment  coefficients  for  the 
bounded  flow  solution  are  then  compared  with  the  free  field 
solution  to  obtain  the  necessary  tunnel  corrections  AC^ ,  and 
ACm  .  The  calculated  drag  coefficient  should  be  zero,  since 
inviscid  flow  is  assumed,  and  serves  as  an  accuracy  check  on 
the  numerical  procedures.  In  the  present  work,  the  free 
stream  values  were  obtained  using  the  singularity  method  of 
Kennedy  [13]  for  the  calculation  of  the  surface  pressure 


yip  ~  Sip 


€ 


d? 


(42) 


A 


. 


'll  ii 
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distribution,  which  was  then  integrated  in  the  same  manner 
as  the  present  numerical  scheme.  Kennedy's  method  was  used 
to  obtain  free  field  values  because  the  present  numerical 
results  were  limited  by  the  maximum  grid  size  which  could  be 
accomodated  in  the  computer  memory,  thus  limiting  the 
distance  to  the  outer  boundary,  and  making’  it  difficult  to 
implement  truly  unbounded  flow  conditions. 

C.  Numerical  Implementation 

The  grid  generation  and  potential  flow  solutions  were 
accomplished  by  iteration  using  successive  over-relaxation 
(SOR).  All  partial  differential  equations  were  discretized 
using  second  order  central  difference  expressions  throughout 
the  solution.  Initial  values  for  the  solutions  were  supplied 
by  linear  interpolation  along  constant  coordinate  lines 
between  known  boundary  values.  Optimum  relaxation  values, 
which  depended  on  several  factors,  were  determined 
empirically  by  observation  of  convergence  rates  as  shown  in 
Figure  21  or  Figure  22.  For  each  grid  type,  optimum 
relaxation  parameters  were  found  for  a  variety  of  grid 
sizes,  geometries,  boundary  point  distributions,  and  grid 
control  amplitude.  Typical  values,  shown  in  Table  3,  are 
slightly  conservative  but  are  applicable  to  a  widerange  of 
grids  of  each  type.  Thompson  [21]  gives  an  algorithm  for 
optimally  varying  the  relaxation  over  the  grid.  This 
algorithm  was  not  used  in  the  present  work,  due  to  the 
capacity  of  the  computer  used,  and  a  constant  relaxation 


' 


85 


Table  3  Relaxation  Parameters  for  Optimum  Convergence 


Grid  Type 


Size  Relaxation  Factor,  W 


0 

61 

100 

X  X 

60 

100 

1.92 

1.94 

c 

100 

X 

30 

1.50 

s 

139 

X 

75 

1.87 

2S 

99 

X 

80 

1.70 
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log(Axm) 


Iteration 


1  argest 
between 


change  on  grid 
i terat i ons 


Figure  21  Effect  of  Relaxation  Parameter  on  Convergence 
Rates  for  an  0-Type  Grid 
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Axloo  =  1  argest  change  on  grid  at  100  iterations 


Figure  22  Determination  of  Optimum  Relaxation  Parameter  for 
an  O-Type  Grid 
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factor  was  used  throughout  the  grid.  Convergence  was  assumed 
when  the  maximimum  change  of  the  solution  on  the  grid  was 
less  than  some  prescribed  value.  Thompson  [24]  showed  that 
for  O-type  grids,  force  coefficients  within  1%  accuracy  were 

-4 

obtained  using  convergence  criteria  of  10  for  grid 
generation  and  1 0 ^  for  potential  flow  solutions.  In  the 
present  work,  these  criteria  were  decreased  an  order  of 
magnitude.  This  did  not  involve  a  substantial  increase  in 
computer  time,  as  seen  in  Figure  21.  However  it  was 
determined  that  for  large  grids  (in  the  physical  space)  the 
single  precision  round-off  error  during  grid  generation  was 
near  5*10“4-  ,  and  the  convergence  criterion  was  adjusted 
accordingly.  Moreover,  the  largest  change  on  the  grid  at 
convergence  occurs  where  the  grid  is  coarsest,  so  the 
relative  error  remained  small. 

The  boundary  conditions  for  the  C-  and  S-type  grids 
lend  themselves  to  solution  by  line  successive 
over-relaxation  (LSOR) .  This  is  a  semi-implicit  method  which 
generates  improved  solution  values  along  an  entire  tj  line 
in  the  computational  grid  by  the  solution  of  a  tri-diagonal 
matrix.  This  is  similar  to  the  alternating  direction 
implicit  method  (ADI),  and  an  attempt  was  made  to  utilize 
the  complete  ADI  method  in  the  solutions.  However,  the 
presence  of  the  cross-derivative  term  in  the  partial 
differential  equations  was  a  complicating  factor:  the 
standard  sequences  of  acceleration  parameters  associated 
with  ADI  solutions  of  elliptic  equations  failed  to  perform 
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as  expected.  Although  convergence  was  initially  accelerated 
in  certain  test  cases,  the  performance  of  the  algorithm  did 
not  justify  the  additional  programming  effort.  Use  of  LSOR 
for  the  same  test  cases  with  the  cross-derivative  terms 
showed  substantial  reductions  in  computer  time  and 
iterations  to  convergence  compared  to  SOR.  Henceforth,  all 
S-type  grid  solutions  were  performed  using  LSOR. 

Some  difficulty  was  encountered  in  the  numerical 
solutions  when  a  high  degree  of  grid  control  was  used.  This 
was  likely  due  to  the  large  error  which  was  sometimes 
present  in  the  initial  solution  values.  During  the  initial 
iterations,  the  grid  adjustments  were  so  large  that  the  grid 
would  fold  back  on  itself  and  the  ensuing  calculations  would 
fail.  It  was  determined  that  the  most  reliable  method  of 
grid  generation  was  then  obtained  when  both  the  relaxation 
parameter  and  grid  control  amplitude  were  initially  set  low 
( W= 1 . 3  and  amplitude  zero)  and  subsequently  increased  in 
steps  as  the  solution  progressed.  These  adjustments  were 
made  automatically  when  the  largest  change  on  the  grid  was 
less  than  some  value,  usually  10  1  .  Once  the  final  values 

for  relaxation  and  grid  control  were  reached  the  solution 
continued  until  convergence  was  obtained. 

It  is  customary  in  numerical  work  to  report  such 
statistics  as  computer  time.  The  University  of  Alberta 
Department  of  Computing  Services  currently  (spring,  1982) 
operates  an  Amdahl  470/V8  mainframe  on  a  time  sharing  basis 
using  the  Michigan  Terminal  System  (MTS).  The  numerical 
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solutions  were  performed  on  a  Floating  Point  Systems  AP-190L 
array  processor,  hosted  by  the  Amdahl.  Exact  timing 
statistics  are  meaningless  unless  accompanied  by  a  lengthy 
discussion  of  the  architecture  and  operating  systems  of  the 
two  machines,  and  the  details  of  the  object  code 
implementation.  As  a  rough  guideline,  the  Amdahl  execution 
speed  is  approximately  6.5  million  instructions  per  second, 
while  the  array  processor,  in  the  configuration  used,  is 
slightly  slower.  Typical  execution  times  were  several 
minutes  on  the  array  processor  and  several  seconds  on  the 
host  Amdahl.  A  large  portion  of  the  host  CPU  time  was  used 
for  disk  I/O  activity.  The  somewhat  slower  array  processor 
was  used  because  of  the  lower  cost  involved,  a  factor  of  40 
less  than  the  Amdahl.  The  array  processor  performed  only  the 
numerical  solutions,  while  all  the  auxiliary  calculations, 
plotting,  and  execution-time  adjustments  were  performed  by 
the  Amdahl.  Some  of  the  O-type  and  C-type  grid  solutions 
could  not  be  implemented  efficiently  due  to  difficulties  in 
operating  the  array  processor,  which  was  under  development 
at  the  time. 

D.  Data  Interpretation 

In  the  words  of  Hamming  [28],  "The  purpose  of  computing 
is  insight,  not  numbers."  In  the  present  work,  where  a 
typical  grid  contained  over  five  thousand  computational 
nodes,  with  up  to  fourteen  variables  defined  at  each  node, 
data  analysis  and  interpretation  was  possible  only  by  means 
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of  the  extensive  graphics  capability  of  the  computing 
system.  Roughly  half  the  host  computer  time  was  used  in  the 
generation  of  graphic  displays  of  the  data.  The  solutions 
sought  here  are  smooth,  due  to  the  nature  of  the  Laplace 
equation,  and  discontinuities  in  any  part  of  the  solution 
became  quickly  evident  when  displayed  graphically.  Suitable 
adjustments  (eg.  the  boundary  point  distribution  and 
implementation  of  more  grid  control)  could  then  be  made  to 
increase  the  overall  suitability  of  the  grid,  and  accuracy 
of  the  solution. 

For  the  grid  generation,  it  was  generally  possible  to 
determine  when  an  adequate  grid  was  obtained  by  merely 
plotting  it,  although  some  very  good  looking  grids  sometimes 
provided  extremely  poor  results.  A  more  quantitative 
approach,  as  suggested  by  Kerlick  and  Klopfer  [22],  involves 
re-arranging  the  metrics  into  the  physically  more  meaningful 
parameters  representing  cell  orientation,  orthogonality, 
cell  aspect  ratio,  and  cell  size.  The  cell  orientation  in 
the  physical  plane,  for  example,  the  rotation  of  the  tangent 
to  an  7)  line  from  the  x  axis,  is  given  by  the  angle  S 
where 

x,  x 

C0S6  =  2  2  1/2  =  1/2 
(x.2  +  y^2)  Y 

The  orthogonality  (or  lack  thereof)  of  the  grid  is 
characterized  by  the  angle  Q  between  the  and  7}  grid 


lines,  where 
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COS@  —  i/o 

(ay)1/2 

The  cell  aspect  ratio,  a,  is  the  ratio  of  the  magnitudes  of 
the  tangent  vectors  to  the  E,  and  77  lines  at  each  grid 
node: 


Finally,  the  Jacobian  represents  the  area  in  the  physical 
space  of  a  single  cell  in  the  computational  space.  One  way 
to  observe  these  properties  is  to  plot  perspective  surface 
drawings  for  the  metrics  and  grid  properties  as  shown  in 
Figures  23  to  28  for  the  S-type  grid  shown  in  Figure  16b.  In 
these  plots  it  is  very  easy  to  observe  areas  of  the  grid 
where  the  properties  vary  rapidly  and  to  make  adjustments  to 
the  grid  accordingly.  Particular  areas  of  importance  are 
near  the  section  leading  and  trailing  edges.  Figures  23  and 
24  indicate  rather  extreme  grid  behavior  at  the  section 
leading  and  trailing  edges.  Figure  26  illustrates  the 
smoothness  of  cell  size  variation  near  the  model.  Figure  27 
indicates  the  maximum  non-orthogonality  (less  than  1  degree 
from  orthogonality)  occurs  well  upstream  of  the  model,  and 
that  near  the  model  the  orthogonality  is  quite  uniform. 
Finally,  Figure  28  shows  that  the  77  -line  spacing  is 
approximately  3  times  the  £,  -line  spacing  near  the  airfoil 
leading  and  trailing  edges. 
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M I N=  0.00122 


MflX=  0.00828 


Figure  23  S-Type  Grid  Alpha  Surface  Plot 


Figure  24  S-Type  Grid  Beta  Surface  Plot 
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Figure  25  S-Type  Grid  Gamma  Surface  Plot 


M I  N=  0.00083 


MRX=  0.01363 


gure  26  S-Type  Grid  Jacobian  Surface  Plot 
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Figure  27 


S-Type  Grid  Orthogonality  Surfao 


Plot 


i 


MIN=  0.19289 
MflX  =  3.57431 


Figure  28  S-Type  Grid  Aspect  Ratio  Surface  Plot 
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Other  graphical  procedures  involved  the  contour 
plotting  of  the  solutions  in  the  physical  plane.  For  the 
stream  function,  these  contours  are  the  flow  streamlines. 
Examples  are  shown  in  Figures  29  to  31  of  the  three  basic 
solutions  for  a  single  element  section,  with  the  total  flow 
shown  in  Figure  32.  Only  part  of  the  solution  field  is  shown 
in  these  Figures:  the  outer  boundary  was  five  chord  lengths 
away.  These  Figures  were  generated  using  an  O-type  grid  and 
the  position  of  the  re-entrant  line  is  indicated,  showing 
the  solution  continuity  across  this  line. 

Other  data  may  be  extracted  from  the  solutions.  Figure 
33  shows  the  pressure  coefficient  contours  through  a  wind 
tunnel  contraction,  constructed  using  two  cubic  profiles  (a 
standard  contraction  shape),  which  clearly  shows  a  velocity 
defect  near  the  tunnel  centre  line  at  the  test  section 
entrance.  Such  data  could  be  useful  in  contraction  design. 
These  examples  show  only  a  sample  of  the  types  of  data  which 
may  be  extracted  from  the  potential  flow  solution  using 
field  techniques-. 

E.  Accuracy  of  BFC  Solutions 

The  O-type  grid  is  most  suited  to  geometries  where  the 
outer  boundary  is  a  large  distance  away  from  the  inner 
boundary.  For  the  wind  tunnel  geometry,  where  the  solution 
domain  is  long  and  thin,  the  O-type  grid  is  very  coarse  and 
highly  skewed  near  the  upstream  and  downstream  boundaries. 
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Figure  29  Stream  Function  Contours  -  Uniform  Stream  at  Zero 
Degrees  (No  Circulation) 


Re-entrant  line 


Figure  30  Stream  Function  Contours  -  Uniform  Stream  at 
Ninety  Degrees  (No  Circulation) 
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Figure  31  Stream  Function  Contours  -  Circulatory  Flow 


Re-entrant  line 


Figure  32  Stream  Function  Contours  -  Total  Flow  at  20 
Degrees  Incidence 
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Figure  33  Pressure  Coefficient  Contours  through  a 
Contraction  Using  a  C-Type  Grid 
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The  C-type  grid,  although  adequate  upstream  of  the  model,  is 
equally  unacceptable  downstream.  Solutions  were  obtained, 
however,  using  the  O-type  grid  for  the  case  of  distant 
boundaries,  to  test  the  validity  of  the  trailing  point  Kutta 
condition,  equation  (38)  and  the  integration  routines  for 
the  calculation  of  the  aerodynamic  coefficients  from  the 
pressure  distribution. 

Thompson  [24]  has  investigated  the  use  of  the  O-type 
grid  for  potential  flow  solutions  for  two  Karman-Tref f tz 
sections.  His  investigation  included  the  effects  of  the 
required  distance  to  "infinity",  the  boundary  point 
distribution,  and  use  of  grid  control,  and  utilized  the 
velocity  Kutta  condition  only.  However,  in  the  present  work, 
it  was  found  impossible  to  implement  the  velocity  condition 
for  a  section  with  a  sharp  trailing  edge,  whether  cusped  or 
wedge  shaped,  using  the  O-type  grid.  This  is  thought  to  be 
due  to  the  extreme  distortion  of  the  grid  at  the  trailing 
edge,  as  shown  in  Figure  34a.  Careful  examination  of  the 
sections  used  by  Thompson  in  [24]  revealed  that  the  sections 
were  modified  to  give  a  rounded  trailing  edge,  by  increasing 
slightly  the  radius  of  the  circle  used  to  generate  the 
sections.  The  grid  around  the  modified  trailing  edge,  shown 
in  Figure  34b,  is  much  more  suitable  for  calculation  of  the 
surface  derivatives  necessary  for  the  velocity  condition 
implementation,  and  explains  the  success  of  Thompson's 
calculations.  Since  the  airfoil  camber  does  not  change  by 
using  the  enlarged  circle,  the  developed  lift  is  the  same. 
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b)  rounded  trailing  edge  of  Thompson ,[ 24 ] 


a)  sharp  trailing  edge 


Figure  34  O-Type  Grid  Trailing  Edge  Details 
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Thompson  indicates  that  results  within  1%  of  the  analytical 
solution  are  readily  obtained. 

The  procedure  of  artificially  rounding  the  trailing 
edge  is  unsuitable  for  the  present  case,  where  the  section 
shape  is  specified.  Calculation  of  the  circulation  using  the 
trailing  Kutta  condition  proved  to  be  a  satisfactory 
alternative  for  these  sharp  trailing  edge  sections  in  an 
O-type  Grid.  Using  a  Joukowski  section,  a  comparison  was 
made  using  the  present  method  with  the  trailing  condition 
and  the  singularity  method  of  Kennedy.  This  was  an  extreme 
case,  as  the  Joukowski  section  was  cusped.  The  same  number 
of  surface  points  were  specified  for  both  calculations, 
results  of  which  are  given  in  Table  4,  along  with  the 
theoretical  values.  It  is  observed  that  the  present  method 
and  the  singularity  method  under-predict  the  lift  at  high 
angle  of  attack.  In  the  range  for  which  unseparated  flow 
might  be  expected  to  occur,  the  lift  coefficients  obtained 
using  the  two  methods  differ  by  less  than  1%,  although  the 
absolute  error  from  the  analytic  solution  is  nearly  4%. 

It  was  difficult  to  assess  the  accuracy  of  the 

solutions  using  the  S-type  grids.  These  results  might  be 

•\ 

expected  to  be  most  accurate  when  the  grid  is  dense  and 
uniform  near  the  model.  Since  the  maximum  grid  size  was 
limited  by  computer  resources,  this  fine  packing  of  the  grid 
was  accomplished  best  by  bringing  the  outer  boundary  towards 
the  model.  This  meant  implementation  of  the  far-field 
boundary  conditions  uncomfortably  close  to  the  model,  and  in 
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Table  4  Comparison  of  Numerical  Solution  with  Singularity 
Solution  of  Kennedy  [13] 


Section:  doukowski  30%  thick,  symmetrical 
30  points  each  surface 

numerical  "infinity"  at  25  chord  lengths 
from  section  mid-chord 


Angle  of 
Attack 
( degrees ) 

Analytic 

solution 

Singulari ty 
solution 
( Kennedy ) 

Numerical 

solution 

0 

o 

o 

o 

o 

0.0012 

2 

0.2699 

0.2603 

0.2616 

4 

0.5395 

0.5203 

0.5217 

6 

0.8085 

0.7797 

0.7811 

8 

1 .0764 

1 .0382 

1 .0396 

10 

1 .3431 

1.2953 

1.2969 

12 

1 .6081 

1.5509 

1.5525 
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fact  represents  the  wind  tunnel  problem  itself,  for  which 
there  are  few  exact  solutions.  The  exact  solution  of 
Tomotika  for  a  flat  plate  in  a  wind  tunnel,  mentioned  by 
Garner  [1],  would  not  include  the  blockage  effect,  which  is 
known  to  interact  with  the  lift  interference.  The  results  of 
de  Vries  and  Schipholt  [15]  could  not  be  checked,  as  the  two 
element  section  profile  used  in  [15]  was  not  specified.  The 
results  of  Williams  and  Parkinson  [2],  which  utilized  a 
singularity  method  for  the  wind  tunnel  problem,  were 
available  for  comparison  with  the  S-type  grid  results.  Using 
a  14%  Clark-Y  airfoil  zero  incidence,  the  lift  was 
calculated  for  a  range  of  tunnel  heights.  The  results, 

Figure  35,  indicate  comparable  accuracy  of  the  finite 
difference  solution  with  the  singularity  method  used  in  [2], 
For  the  wind  tunnel  problem,  solution  accuracy  was 
dependent  on  the  extent  of  grid  control.  The  O-type 
solutions  were  much  more  sensitive  in  this  respect  than  the 
S-type  grids  and  for  the  tunnel  geometries  the  S-type  grids 
were  used  exclusively.  With  experience,  the  degree  of  grid 
control  necessary  for  an  accurate  solution,  which  depended 
on  several  parameters  (grid  size,  tunnel  geometry,  section 
profile,  etc.),  was  determined  by  trial  and  error,  using  the 
metric  plots  described  previously. 
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Figure  35  Comparison  of  S-Type  Solution  With  Singularity 
Solution  for  a  Wind  Tunnel  Geometry 


F.  Conclusion 


For  the  present  application  to  the  calculation  of  wind 

tunnel  corrections,  the  overall  method  is  quite  extravagent 

in  terms  of  computational  effort.  One  of  the  purposes  of  the 

present  work  is  to  assess  the  accuracy  of  the  method  of 
,  » 

boundary  fitted  coordinates.  Application  to  the  real  problem 
of  determining  the  high  blockage  tunnel  corrections  provides 
a  suitable  and  practical  test  case.  Obviously,  as  with  most 
extensive  computational  efforts,  the  immediate  application 
of  the  algorithm  by  other  users  cannot  be  expected  (a 
manifestation  of  the  "not  invented  here"  syndrome).  The 
computational  program  is  quite  large  and  expensive  to  run. 
Extensive  use  of  the  present  procedures  would  result  in  a 
series  of  charts,  tables,  or  correction  relations  which 
would  be  of  more  immediate  practical  use  at  other  test 
installations.  This  data  base  would  be  quite  extensive, 
because  of  the  number  of  parameters  involved.  The  generation 
of  this  extended  data  base  was  not  a  feasible  undertaking 
due  to  time  and  fiscal  constraints. 


IV.  Experiments 


To  investigate  the  applicability  of  the  proposed  wall 
interference  corrections  calculated  using  the  methods  of 
Chapter  III,  experiments  were  carried  out  in  the  University 
of  Alberta  low  speed  wind  tunnel  on  several  models  using 
different  tunnel  geometries.  All  tests  were  performed  at  a 
chord  Reynolds  number  of  106  .  Two  single  element  airfoils 
and  one  two  element  airfoil  were  tested.  The  tunnel  height 
could  be  adjusted  by  the  movement  of  two  parallel  interior 
wall  inserts,  one  on  each  side  of  the  two  dimensional 
models.  Pressure  data  was  acquired  using  several  high 
sensitivity  differential  pressure  transducers  and  a  computer 
controlled  data  acquisition  system.  Pressure  distributions 
were  integrated  using  the  trapezoidal  rule  to  give  the 
uncorrected  force  and  moment  coefficients. 

Free  field  experimental  data  for  all  three  models  was 
obtained  using  the  standard  low  blockage  corrections 
described  in  Chapter  II  on  data  obtained  with  the  inserts 
removed,  giving  a  chord-to-height  ratio  of  0.41.  For  these 
tests,  the  section  drag  coefficient  was  obtained  using  a 
wake  rake. 

The  interference  tests  involved  several  wall  positions 
yielding  a  broad  range  of  chord-to-height  ratios.  The  tunnel 
wall  pressure  distribution  was  measured  for  comparison  with 
the  inviscid  theory.  Due  to  physical  constraints,  no 
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provision  was  made  for  drag  measurements  during  the  blockage 
tests . 

Data  for  each  tunnel  configuration  was  obtained  for  a 
range  of  angles  of  attack  where  attached  flow  on  the 
sections  was  maintained  since  the  present  potential  flow 
theory  is  based  on  the  assumption  that  viscous  effects  are 
neglected.  The  final  data  for  the  interference  tests  was 
then  reduced  using  both  the  present  corrections  and  the 
standard  constant  incidence  corrections  for  comparison  with 
the  free  field  tests. 

A.  Test  Equipment,  Procedure,  and  Data  Reduction 

Tunnel  Calibration 

The  University  of  Alberta  Department  of  Mechanical 
Engineering  operates  a  closed  circuit,  single  return,  low 
speed  wind  tunnel  with  a  test  section  1.2m  high,  2.4m  wide 
and  11.0m  in  length.  The  tunnel  can  be  operated  as  a  low 
speed  boundary  layer  tunnel  or  high  speed  aerodynamic 
tunnel,  with  precise  speed  control.  For  aerodynamic  testing, 
the  maximum  unit  Reynolds  number  is  2.0*106  per  meter, 
depending  on  the  tunnel  temperature,  atmospheric  pressure, 
and  the  presence  of  filters  upstream  of  the  settling 
chamber.  A  turbulence  level  in  the  test  section  of  less  than 
0.1%  is  obtained  using  three  fine  mesh  screens  at  the 
entrance  to  the  6.2:1  contraction.  Two  dimensional  models, 
usually  of  one  meter  chord,  are  mounted  vertically  in  the 
test  section,  having  therefore  a  span  of  1.2  meters. 
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The  tunnel  airspeed  was  determined  using  the  procedure 
described  in  Pope  and  Harper  [5].  The  change  in  static 
pressure  between  the  settling  chamber  and  a  location 
approximately  three-quarters  of  the  length  through  the 
contraction  was  measured.  This  was  calibrated  against  the 
test  section  dynamic  pressure  at  the  model  location  in  the 
empty  tunnel  obtained  using  a  standard  pitot-static  probe, 
corrected  for  stem  and  tip  error.  This  procedure  was  adopted 
when  preliminary  tests  indicated  that  the  static  pressure 
field  generated  by  the  large  models  extended  quite  far 
upstream,  making  local  measurement  of  the  tunnel  free  stream 
static  pressure  unreliable.  The  final  airspeed  calibration 
was  insensitive  to  model  orientation. 

For  the  blockage  tests,  the  test  section  dynamic 
pressure  was  determined  by  a  pitot-static  probe  mounted 
between  the  wall  inserts,  approximately  two  chord  lengths 
upstream  from  the  model  leading  edge.  Tunnel  speed  was  set 
for  a  low  angle  of  attack  for  which  the  upstream  influence 
of  the  model  is  small,  after  which  the  probe  was  removed 
from  the  flow.  The  precise  speed  control  of  the  tunnel  drive 
was  thereafter  able  to  maintain  the  tunnel  speed  within  1%, 
as  monitored  by  the  differential  pressure  through  the 
contraction,  provided  that  flow  separation  on  the  model  was 
avoided,  which  would  increase  the  total  tunnel  loss 
coefficient  and  retard  the  flow.  This  procedure  was 
necessary  since  the  models  used  were  extremely  sensitive  to 
wake  turbulence  produced  by  the  probe.  Subsequently,  the 
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test  section  free  stream  static  pressure  was  determined  from 
measurements  near  the  leading  edge  of  the  wall  inserts, 
where  the  indicated  pressure  was  observed  to  be  the  same  on 
both  walls. 


Models 

To  ensure  the  coverage  of  a  wide  range  of  test 
geometries,  three  airfoil  sections  were  tested.  Section 
profiles  for  the  three  models  are  shown  in  Figure  36. 

The  first  model  was  a  single  element  high  performance 
section  which  had  been  designed  for  very  low  drag  (Cd=.007) 
at  the  test  Reynolds  number.  This  minimized  the  wake 
blockage  and  boundary  layer  displacement  effects,  allowing 
closer  matching  of  the  real  flow  with  the  theoretical  flow, 
so  long  as  separation  was  avoided. 

The  second  model,  also  a  single  element  section,  was 
chosen  for  its  very  thick  cross-section  and  high  camber. 

This  section,  designed  by  Kennedy  [13],  develops  the  highest 
lift  coefficient  for  a  single  element  section  known  to  the 
present  author.  Unfortunately,  as  opposed  to  the  first 
model,  the  viscous  effects  are  very  large  and  the  boundary 
layer  and  wake  effects  are  not  negligible.  At  a  Reynolds 
number  of  106  ,  Kennedy  found  that  the  combined  displacement 
thicknesses  of  the  upper  and  lower  surface  boundary  layers 
at  the  trailing  edge  was  23mm  for  the  design  angle  of  15 
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degrees.  The  drag  coefficient  was  correspondingly  high, 
Cd=0.0284,  and  the  maximum  lift  obtained  was  significantly 
less  than  the  potential  flow  value.  Nevertheless,  the 
lift/drag  ratio  was  still  quite  high,  and  the  lift 
interference  and  solid  blockage  corrections  predominated. 

The  third  section  tested  was  a  two  element  section 
based  on  model  #1  with  a  35%  chord  slotted  flap,  deflected 
20  degrees.  This  section  develops  quite  high  lift/drag 
ratios,  and  again  allows  the  real  flow  to  correspond  closely 
to  the  inviscid  theory. 

The  boundary  layer  structure  on  each  of  these  three 
models  was  relatively  insensitive  to  changes  in  the  pressure 
distributions  caused  by  incidence  changes  when  operated 
within  the  laminar  drag  bucket.  Pope  and  Harper  [5]  suggest 
that  changes  in  boundary  layer  structure  caused  by 
distortion  of  the  pressure  distribution  becomes  important 
only  for  very  high  c/h  ratios.  It  was  expected,  then,  that 
distortion  of  the  pressure  distributions  during  the  blockage 
tests  would  not  produce  large  changes  in  the  section 
performance.  This  was  subsequently  confirmed  for  moderate 
blockage  values.  For  the  tests  at  the  highest  blockage 
values,  the  effect  of  the  wall-induced  pressure  distribution 
distortion  became  quite  pronounced.  Each  of  the  single 
element  sections  utilizes  a  Stratford-type  turbulent 
boundary  layer  pressure  recovery,  which  allows  maximum 
pressure  recovery  with  minimum  drag.  The  Stratford  boundary 
layer  develops  with  a  constant  shape  factor  just  on  the 
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Figure  36  Profiles  of  Models  Tested 


point  of  separation.  At  high • blockage ,  the  distortion  of  the 
pressure  distribution  disrupted  this  recovery,  resulting  in 
changes  in  the  section  performance.  Particularly,  for  model 
#2,  the  maximum  lift  was  affected,  as  well  as  performance  of 
the  section  leading  up  to  stall. 

i 

All  models/were  fitted  with  0.5mm  diameter  static 
pressure  taps  around  the  section  on  the  model  centre-line, 
concentrated  in  regions  of  high  pressure  gradients.  Static 
pressure  tubes  were  led  out  of  the  tunnel  through  the  model 
interior  and  connected  to  a  multi-tube  scanivalve  which 
sequentially  connected  each  tube  to  a  sensitive  pressure 
transducer . 

The  models  were  mounted  on  3/16"  thick  aluminum 
endplates,  equipped  with  boundary  layer  suction  holes 
adjacent  to  the  model  surface  in  regions  of  adverse  pressure 
gradients.  The  model  and  endplates  were  then  mounted  on 
circular  turntables  in  the  tunnel  ceiling  and  floor  which 
could  be  rotated  to  give  the  desired  angle  of  attack.  Each 
turntable  was  fitted  with  externally  mounted  suction  plenum 
chambers,  which  were  connected  to  a  large  centrifugal  fan 
through  6"  diameter  rigid  plastic  pipe.  The  plenum  pressure 
was  maintained  at  2.8kPa  below  atmospheric  pressure 
(Cp=-10).  The  purpose  of  the  boundary  layer  suction  was  to 
remove  portions  of  the  thick  boundary  layer  on  the  wind 
tunnel  floor  and  ceiling.  In  the  presence  of  adverse 
pressure  gradients  generated  by  the  model,  separation  of 
these  boundary  layers  would  seriously  degrade  the  two 
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dimensional  character  of  the  flow  by  allowing  formation  of 
strong  trailing  vortices. 


Wall  Inserts 

The  wall  inserts  were  made  up  of  3/4"  plywood  sheets. 
Initial  tests  were  carried  out  on  the  low  drag  section  (#1) 
with  the  inserts  extending  2.44  chord  lengths  fore  and  aft 
of  the  model  mid-chord.  This  was  subsequently  increased  to 
3.05  chord  lengths,  for  a  total  extent  of  the  inserts  of 
6.10  meters.  The  inserts  extended  forward  to  the  tunnel 
contraction  exit  plane.  Each  insert  panel  was  fitted  with  24 
1/8"  diameter  static  pressure  taps  at  a  spacing  of  6"  near 
the  model  and  12"  at  the  extremities.  The  total  extent  of 
the  wall  pressure  tap  coverage  was  4.42  chord  lengths.  Each 
tap  was  carefully  flush  mounted  and  counter-sunk  to  ensure 
accurate  static  pressure  measurements.  The  inserts  were 
attached  to  the  floor  and  ceiling  of  the  tunnel,  and  were 
kept  parallel  within  +/-  3mm  over  the  entire  length.  At  the 
low  dynamic  pressure  of  the  tests,  deflection  of  the  inserts 
was  negligible.  The  leading  edge  of  the  inserts  was  rounded 
to  ensure  smooth  entry  of  the  flow  into  the  reduced  test 
section.  During  preliminary  runs,  vertical  wires  with 
numerous  6"  tufts  were  placed  upstream  of  each  leading  edge 
to  observe  any  unwanted  flow  curvature.  No  serious  flow 
irregularities  were  observed  and  the  wires  were  later 
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removed . 

The  test  geometry  for  the  free  field  tests  is  shown  in 
Figure  37.  The  wake  rake  apparatus  is  shown  downstream  of 
the  model.  Figure  38  shows  the  same  model  with  the  wall 
inserts  in  place  at  a  chord  to  height  ratio  of  0.93,  the 
maximum  possible  due  to  overlapping  of  .the  inserts  with  the 
model  turntables. 

Pressure  Measurements  -  Transducers  and  Calibrations 

Differential  pressures  were  measured  using  several 
Validyne  pressure  transducers,  all  supplied  with  the  tunnel 
stagnation  pressure  as  reference.  The  transducers  were 
equipped  with  appropriate  diaphragms  (1"  water,  0.1  psi  or 
0.2  psi)  and  calibrated  against  an  accurate  inclined 
manometer.  All  calibrations  were  linear,  and  were  re-checked 
approximately  every  10  hours  of  running  time.  The  maximum 
drift  observed  in  the  sensitivity  of  the  transducers  was 
less  than  0.5%.  The  data  acquisition  system  automatically 
checked  and  adjusted  the  calibration  zeros  at  the  start  of 
the  data  acquisition  sequence,  when  all  transducers  were 
supplied  with  zero  pressure  differential.  To  diminish  the 
effect  of  random  pressure  fluctuations,  each  measurement  was 
repeated  three  to  five  times  and  averaged.  In  addition,  when 
using  the  wake  rake  for  drag  measurements,  the  wake  pressure 
transducer  output  signal  was  sent  through  a  low  pass  filter 
(cut  off  0.3  hz . ) ,  to  smooth  out  fluctuations  due  to  wake 
turbulence.  During  the  data  sequence,  the  scanivalve 
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Figure  37  Tunnel  Geometry  for  Free  Field  Tests 
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Figure  38  Tunnel  Geometry  for  Blockage  Tests 
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pressure  transducers  were  allowed  to  stabilize  for  1/2 
second  before  readings  were  taken. 

Data  Acquisition  and  Reduction 

The  large  extent  of  these  experimental  investigations 
would  not  have  been  possible  without  automated  data 
acquisition  and  reduction.  The  system  configuration  is  shown 
in  Figure  39.  During  the  acquisition  phase,  under  control  of 
the  central  computer  (Hewlett-Packard  HP-85),  the  scanner 
would  route  the  desired  pressure  transducer  signal  to  the 
digital  voltmeter.  Voltage  readings  were  returned  to  the 
computer  and  converted  to  pressures  using  the  appropriate 
calibration.  The  scanivalves  would  then  be  stepped  to  the 
next  location  and  measurements  repeated  until  all  pressure 
taps  had  been  monitored.  The  resulting  pressure 
distributions  were  then  integrated  using  the  trapezoidal 
rule  to  obtain  the  model  force  and  moment  coefficients. 
Finally,  the  pressure  distributions  could  be  plotted,  either 
in  soft  copy  on  the  computer  screen,  or  in  hard  copy  on  the 
plotter.  All  pressure  data  was  stored  on  magnetic  tape 
cartridges  for  later  retrieval  and  further  analysis  on  the 
University  main  computer,  if  desired.  The  entire  data 
acquisition/reduction  sequence  for  a  single  angle  of  attack 
could  be  accomplished  in  under  three  minutes  with  all 
options  in  effect.  For  the  free  field  tests,  the  standard 
corrections  were  incorporated  in  the  system  program.  At  the 
conclusion  of  a  run,  a  lift  curve  (and  drag  polar  for  free 
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Transducer  signals 


Figure  39  Data  Acquisition  Schematic 
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field  tests)  was  automatically  produced. 

Free  Field  Drag  Calculation 

One  of  the  most  difficult  quantities  to  measure  is  the 
section  profile  drag.  The  procedure  adopted  here  is  based  on 
the  Jones  formula  (Schlicting  [29])  for  the  determination  of 
the  profile  drag  when  the  static  pressure  at  the  wake 
measurement  location  is  unequal  to  the  free  stream  static 
pressure.  It  was  assumed  that  the  stagnation  pressure 
outside  the  wake  was  the  same  as  the  stagnation  pressure 
ahead  of  the  model,  and  that  the  static  pressure  was 
constant  through  the  wake.  The  resulting  formula  is 


(44) 


where  q-j  is  the  main  tunnel  dynamic  pressure,  calculated 
from  the  tunnel  contraction,  is  the  dynamic  pressure 

outside  the  wake,  and  Ap  is  the  drop  in  stagnation  pressure 
within  the  wake  as  measured  by  the  transducer.  The  integral 
was  evaluated  using  the  trapezoidal  rule  across  the  width  of 
the  wake. 

B.  Free  Field  Results 

The  raw  and  corrected  experimental  results  for  the 
three  models  are  shown  in  Figures  40  to  42  for  the  lowest 
c/h  ratio  of  the  tests  (c/h=0.41).  Both  the  first  and  second 
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Figure  40  Model  #1 


Free  Field  Lift  and  Moment  Coefficients 
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Figure  41  Model  #2  -  Free  Field  Lift  and  Moment  Coefficients 


o  o 


.50  1.00  1.50  2.00 


128 


o 

o 


„  C 


LO 

o 

i 


0 


Figure  42  Model  #3 


Free  Field  Lift  and  Moment  Coefficients 


129 


order  corrections  were  made,  using  the  constant  incidence 
forms,  equations  (20)  and  (21).  The  lift  curve  slope  used  in 
the  constant  incidence  corrections  was  evaluated  using  the 
corrected  data  in  each  case  using  the  iterative  procedure 
described  previously.  Both  corrections  yield  similar 
results,  with  the  first  order  lift  corrections  being 
slightly  larger  than  the  second  order  forms.  The  two  forms 
for  the  moment  corrections  yield  practically  identical 
results.  For  this  c/h  ratio,  the  first  order  corrections  are 
probably  inappropriate  and  it  was  assumed,  therefore,  that 
the  results  obtained  using  the  second  order  forms  represent 
true  free  field  values. 

Figures  43  to  45  show  the  contributions  of  wake 
blockage,  solid  blockage,  and  lift  interference  to  the  total 
lift  correction  for  each  model  for  these  free  field  results. 
The  wake  blockage  correction  is  a  minor  portion  of  the  total 
(typically  less  than  10%)  and  for  the  wall  interference 
tests  was  therefore  ignored. 


C.  High  Blockage  Results 

The  effect  of  the  wall  position  on  the  uncorrected  data 
is  shown  in  Figures  46  to  48  for  the  three  models.  The  free 
field  results  are  also  included  in  these  figures  and  it  is 
observed  that  at  high  blockage  and  lift  coefficient,  the 
lift  correction  approaches  30%  of  the  experimental  value. 
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Figure  43  Contributions  of  Corrections  -  Model  #1,  Free 
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Figure  44  Contributions  of  Corrections  -  Model  #2,  Free 
Field 
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Figure  45  Contributions  of  Corrections  -  Model  # 3,  Free 
Field 
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Figure  46  Model  #1  -  Effect  of  Wall  Position  on  Uncorrected 


Data 
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Figure  47  Model  #2  -  Effect  of  Wall  Position  on  Uncorrected 


Data 
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Figure  48  Model  #3  -  Effect  of  Wall  Position  on  Uncorrected 


Data 
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This  high  blockage  data  was  corrected  using  the 
standard  forms,  but  ignoring  the  wake  blockage,  to  yield  the 
results  shown  in  Figures  49  to  51.  The  first  order 
corrections  are  consistently  too  large,  while  the  second 
order  corrections  are  too  small,  and  the  errors  in  these 
corrections  increase  with  tunnel  blockage  and  lift.  The 
standard  forms,  both  first  and  second  order,  are  clearly 
unsatisfactory  for  c/h  ratios  greater  than  0.5. 

The  present  corrections  were  implemented  as  follows. 

The  theoretical  analysis  was  used  to  generate  correction 
tables  for  each  model  which  contained  values  of  the  lift  and 
moment  increments  as  functions  of  the  uncorrected  lift 
coefficient  and  c/h  ratio.  The  correction  tables  for  the 
three  models  are  displayed  graphically  in  Figures  52  to  54. 
For  a  given  test,  the  relevant  correction  was  obtained  by 
entering  the  table  at  the  experimentally  obtained  lift 
coefficient  and  c/h  ratio  and  extracting  the  relevant 
correction.  This  procedure  was  automated  using  a  two-way 
cubic  spline  interpolation,  shown  schematically  in  Figure 

55.  The  correction  table  data  was  used  to  construct  the 
splines  SI  to  S5  for  five  c/h  ratios  between  0.2  and  1.0, 
which  bracketed  the  experimental  tests.  The  increment  at  the 
uncorrected  lift  coefficient  was  then  determined  by 
interpolation  on  SI  to  S5.  A  cross-spline  was  then  fitted, 

56,  to  these  increments  and  the  final  correction  was 
extracted  by  interpolation  on  S6  at  the  experimental  c/h 
ratio.  This  somewhat  involved  procedure  allowed  a  great 
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Figure  49  Model  #1  -  High  Blockage  Results  (Standard 
Corrections) 
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Figure  50  Model  #2  -  High  Blockage  Results  (Standard 
Corrections ) 
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Figure  51  Model  #3  -  High  Blockage  Results  (Standard 
Corrections ) 
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Figure  52  Model  #1  -  Correction  Table 


141 


Figure  53  Model  #2  -  Correction  Table 


Figure  54  Model  #3  -  Correction  Table 
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Figure  55  Spline  Interpolation  Scheme  for  Correction  Tables 
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reduction  in  the  number  of  complete  theoretical • calculations 
used  to  generate  the  table:  the  data  base  consisted  of 
calculations  at  five  different  angles  of  attack,  for  each  of 
the  five  c/h  ratios.  The  splines  SI  to  S5  were  evaluated 
only  once  for  each  model. 

% 

The  data  corrected  using  the  above  procedure  is  shown 
in  Figures  56  to  58  for  the  three  models,  along  with  the 
free  field  results  obained  using  the  second  order 
corrections  for  the  lowest  c/h  ratio.  There  is  good 
agreement  of  the  corrected  data  with  the  free  field  data. 
Results  for  model  #1,  a  more  typical  airfoil  shape,  are 
excellent.  The  results  for  model  #2  are  satisfactory  for  all 
but  the  highest  c/h  ratio,  where  the  effect  of  the  pressure 
distribution  distortion  is  evident.  The  two-element  section 
lift  results  are  satisfactory;  however  the  moment 
corrections  are  not  an  improvement  over  the  standard  forms. 
This  was  thought  to  be  due  to  errors  in  the  theoretical 
pressure  distribution  in  the  slot  region  on  the  main 
section,  caused  by  the  grid  geometry.  Although  not  largely 
affecting  the  lift,  this  had  a  large  effect  in  the 
theoretical  pitching  moment. 

For  a  normal  test,  at  a  single  c/h  ratio,  the  above 
procedure  using  the  complete  correction  table  is  not 
required.  In  this  case,  theoretical  lift  and  moment 
increments  may  be  calculated  for  several  angles  of  attack 
for  the  experimental  c/h  ratio.  A  single  spline  may  then  be 
fitted  through  this  data,  and  the  correction  interpolated  at 
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Figure  56  Model  #1  -  High  Blockage  Results  Using  Correction 
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Figure  57  Model  #2  -  High  Blockage  Results  Using  Correction 
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Figure  58  Model  #3  -  High  Blockage  Results  Using  Correction 
Table 
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the  experimental  lift  coefficient. 


D.  Conclusion 

The  numerically  generated  corrections  are  a 
satisfactory  improvement  over  the  standard  forms  for  the 
three  models  tested.  This  indicates  the  possibility  of 
utilizing  a  wider  range  of  tunnel  geometries,  in  particular 
higher  c/h  ratios,  for  two-dimensional  testing.  The  effect 
of  pressure  distribution  distortion  which  occurs  at  very 
high  blockage  places  an  upper  limit  on  the  allowable  range 
of  c/h  ratios  for  which  reliable  corrections  can  be  made. 
This  upper  limit  is  near  c/h=0.7,  nearly  twice  the  range  for 
which  the  standard  corrections  are  considered  valid.  The 
numerical  corrections  apply  equally  well  for  the  two-element 
section,  although  great  care  in  obtaining  the  theoretical 
solution  is  necessary.  These  results  indicate  that  the 
present  correction  technique  is  valid,  and  can  be  used  for  a 
wide  variety  of  wing  profiles  and  wind  tunnel  geometries. 


V.  Conclusion 


Aerodynamic  performance  data  obtained  in  wind  tunnels 
is  subject  to  the  influence  of  the  tunnel  boundaries.  To 
obtain  free  field  data,  the  experimental  results  must  be 
corrected  for  this  wall  interference.  For  models  which  are 
small  in  relation  to  the  tunnel  height,  the  standard 
correction  procedures  are  adequate;  however,  if  the  model  is 
appreciably  larger,  the  standard  corrections  are  no  longer 
suitable.  A  method  has  been  developed  which  will  accurately 
predict  the  lift  and  pitching  moment  corrections  for  these 
extreme  geometries. 

This  estimation  of  the  wall  interference  is  based  on  an 
accurate  potential  flow  solution  for  the  flow  in  the  tunnel. 
The  effect  of  the  walls  is  contained  explicitly  in  the 
analysis.  There  are  no  restrictions  on  the  model  size  or 
tunnel  geometry,  and  the  method  applies  equally  well  to 
multi-element  sections.  The  solution  contains  the  coupling 
of  the  lift  interference  and  solid  blockage  effects,  which 
is  not  accounted  for  in  the  standard  correction  procedures. 

The  method  utilizes  numerically  generated  boundary 
fitted  coordinate  systems  which  allow  accurate  and  efficient 
numerical  solutions  for  regions  with  irregularly  shaped 
boundaries.  The  coordinate  system  is  obtained  by  the  finite 
difference  solution  of  a  system  of  elliptic  partial 
differential  equations,  using  either  successive 
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over-relaxation  (SOR)  or  line  successive  over-relaxation 
(LSOR)  where  appropriate.  Several  grid  topologies  have  been 
investigated  for  application  to  the  problem  of  determining 
the  wind  tunnel  corrections  for  two-dimensional 
incompressible  flow.  This  investigation  resulted  in 
determination  of  optimal  numerical  procedures  for  each  grid 
type.  Use  of  LSOR,  where  possible,  was  extremely  beneficial 
in  reducing  computation  time  and  cost. 

Some  care  is  required  in  choosing  a  grid  which  will 
lead  to  accurate  results.  The  most  appropriate  grid  geometry 
for  the  wind  tunnel  problem  is  the  slit  type,  due  to  the 
uniformity  of  the  grid  properties  throughout  the  solution 
field.  For  the  multi-element  section,  using  the  double-slit 
geometry,  adequate  results  can  be  obtained  only  by  extensive 
use  of  special  grid  control  functions.  The  decaying 
exponential  form  of  the  control  function  is  easily 
implemented  and  quite  adequate  for  the  cases  treated  here. 

The  potential  flow  solution  for  the  stream  function  is 
obtained  using  similar  numerical  procedures  as  the  grid 
generation.  The  lifting  potential  flow  is  found  by 
superposition  of  the  basic  solutions  for  uniform  onset  flow 
and  pure  circulatory  flow.  The  circulation  of  the  lifting 
flow  is  set  by  the  Kutta  condition.  Implementation  of  the 
Kutta  condition  in  the  numerical  scheme  was  investigated  in 
detail:  both  trailing  point  and  trailing  edge  velocity 
conditions  were  examined  in  two  basic  grid  types.  For  the  0- 
and  C-type  grids,  the  trailing  point  condition  was  superior 
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due  to  ease  of  implementation  and  the  ability  to  treat  sharp 
trailing  edges.  For  the  slit  geometries,  the  velocity 
condition  was  superior. 

The  potential  flow  numerical  solution  compares  in 
accuracy  with  other  analysis  techniques  when  suitable  grids 
are  used.  Comparison  of  the  numerical  results  for  the  wind 
tunnel  geometry  with  free  field  results  obtained  using  a 
proven  singularity  method  yielded  the  desired  wall 
corrections. 

To  test  the  accuracy  of  the  numerically  generated 
corrections,  wind  tunnel  tests  were  carried  out  on  three 
two-dimensional  wing  models,  including  one  two-element 
section,  of  markedly  different  geometry  and  performance. 
Movement  of  two  large  wall  inserts  within  the  tunnel  allowed 
testing  for  a  wide  range  of  tunnel  geometries.  Wind  tunnel 
data  was  obtained  using  an  accurate  and  efficient  data 
acquisition  system. 

The  measured  wind  tunnel  data  was  corrected  using  both 
the  standard  correction  forms  and  the  numerically  generated 
corrections.  The  latter  proved  to  be  far  superior  to  the 
standard  forms  for  extreme  tunnel  geometries.  For  all  three 
models  the  numerical  corrections  were  able  to  reduce  the 
measured  data  for  a  wide  range  of  tunnel  geometries  to  a 
common  set  of  data  indicative  of  the  free  field  performance. 
The  corrections  are  most  notably  accurate  for  the  single 
element  high  lift  section  and  the  two-element  section.  For 
very  extreme  geometries,  however,  the  corrections  become 


' 


152 


somewhat  in  error  due  to  the  effect  of  the  tunnel  walls  on 
the  wing  boundary  layers  through  distortion  of  the  airfoil 
pressure  distribution.  In  this  case,  it  is  doubtful  that  any 
potential  flow  based  wall  correction  scheme  is  suitable. 
Models  with  such  sensitive  boundary  layer  structure  impose 
special  limitations  on  the  allowable  test  geometries.  In 
general,  accurate  corrections  may  be  obtained  for  model 
chord  to  tunnel  height  ratios  up  to  0.7,  nearly  double  the 
range  of  the  standard  correction  forms. 

The  numerical  technique  using  boundary  fitted 
coordinates  is  expensive  and  requires  a  substantial  computer 
capacity.  In  addition,  some  subjective  input  is  required  to 
generate  accurate  solutions.  In  this  regard,  the  need  for 
automatic  control  of  the  grid  properties  (the  adaptive 
grids)  is  well  demonstrated.  The  final  solution,  however, 
contains  a  wealth  of  easily  accessible  information.  For  the 
wind  tunnel  problem,  the  majority  of  this  information  is  not 
used.  Nonetheless,  the  results  obtained  indicate  that 
acceptable  accuracy  of  the  solution  exists  and  that  finite 
difference  potential  flow  solutions  can  be  an  accurate  and 
useful  analysis  tool. 
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